#This is the replication file for Thomann, E. van Engen, N. and L. Tummers. Forthcoming. The necessity of discretion: a behavioral evaluation of bottom-up implementation theory.  Journal of Public Administration Research and Theory.

# This code replicates the analysis of dataset 2, direct calibration method.

#R code based  on Thomann, Eva and Stefan Wittwer (2016). Performing fuzzy- and crisp set QCA with R: A user-oriented beginner's guide. Version 19.9.2016. URL: http://www.evathomann.com/links/qca-r-manual [date of access: 25.11.2016].


#clear workspace
rm(list = ls())


library(lattice); library(arm); library(xtable); library(foreign); library(psych); library(directlabels); library(betareg); library(VIM); library(base); library(plyr); library(dplyr); library(QCA);  library(SetMethods)

setwd("C:/Users/Eva/Dropbox/Arbeit/Publikationen/Papers/Two-factor theory/Material/Data/Dataset 2")


##############PREPARING THE DATASET###########

tfr <- read.csv("tfr_d2.csv", header = TRUE, sep = ";",  dec = ",")
describe(tfr)

#convert -99 into missings
tfr[tfr==-99] <- NA
describe(tfr)

#check missings

colnames(tfr)
describe(tfr)

missings <- as.numeric(complete.cases(tfr))
tabmiss <- table(missings)
prop.table(tabmiss)

#identify  cases with missings on all indicators for a causal factor
##for W
colnames(tfr)

tfr$w1miss <- as.numeric(is.na(tfr$W1))
tabw1miss <- table(tfr$w1miss)
round(prop.table(tabw1miss), digits=3)


tfr$w2miss <- as.numeric(is.na(tfr$W2))
tabw2miss <- table(tfr$w2miss)
round(prop.table(tabw2miss), digits=3)

tfr$w3miss <- as.numeric(is.na(tfr$W3))
tabw3miss <- table(tfr$w3miss)
round(prop.table(tabw3miss), digits=3)

tfr$w4miss <- as.numeric(is.na(tfr$W4))
tabw4miss <- table(tfr$w4miss)
round(prop.table(tabw4miss), digits=3)

tfr$w5miss <- as.numeric(is.na(tfr$W5))
tabw5miss <- table(tfr$w5miss)
round(prop.table(tabw5miss), digits=3)

wmiss <- subset(tfr, select = c("w1miss", "w2miss", "w3miss", "w4miss", "w5miss"))
tfr$wmiss <- rowSums(wmiss)
tfr$wmiss
tfr$wmissout <- as.numeric(tfr$wmiss ==5)
tfr$wmissout
tabwmissout <- table(tfr$wmissout)
prop.table(tabwmissout)


##for SP

tfr$sp1miss <- as.numeric(is.na(tfr$SP1))
tabsp1miss <- table(tfr$sp1miss)
round(prop.table(tabsp1miss), digits=3)


tfr$sp2miss <- as.numeric(is.na(tfr$SP2))
tabsp2miss <- table(tfr$sp2miss)
round(prop.table(tabsp2miss), digits=3)

tfr$sp3miss <- as.numeric(is.na(tfr$SP3))
tabsp3miss <- table(tfr$sp3miss)
round(prop.table(tabsp3miss), digits=3)

tfr$sp4miss <- as.numeric(is.na(tfr$SP4))
tabsp4miss <- table(tfr$sp4miss)
round(prop.table(tabsp4miss), digits=3)

tfr$sp5miss <- as.numeric(is.na(tfr$SP5))
tabsp5miss <- table(tfr$sp5miss)
round(prop.table(tabsp5miss), digits=3)

tfr$sp6miss <- as.numeric(is.na(tfr$SP6))
tabsp6miss <- table(tfr$sp6miss)
round(prop.table(tabsp6miss), digits=3)

spmiss <- subset(tfr, select = c("sp1miss", "sp2miss", "sp3miss", "sp4miss", "sp5miss", "sp6miss"))
tfr$spmiss <- rowSums(spmiss)
tfr$spmiss
tfr$spmissout <- as.numeric(tfr$spmiss==6)
tfr$spmissout
tabspmissout <- table(tfr$spmissout)
prop.table(tabspmissout)

##for TP
tfr$tp1miss <- as.numeric(is.na(tfr$TP1))
tabtp1miss <- table(tfr$tp1miss)
round(prop.table(tabtp1miss), digits=3)


tfr$tp2miss <- as.numeric(is.na(tfr$TP2))
tabtp2miss <- table(tfr$tp2miss)
round(prop.table(tabtp2miss), digits=3)

tfr$tp3miss <- as.numeric(is.na(tfr$TP3))
tabtp3miss <- table(tfr$tp3miss)
round(prop.table(tabtp3miss), digits=3)

tfr$tp4miss <- as.numeric(is.na(tfr$TP4))
tabtp4miss <- table(tfr$tp4miss)
round(prop.table(tabtp4miss), digits=3)

tfr$tp5miss <- as.numeric(is.na(tfr$TP5))
tabtp5miss <- table(tfr$tp5miss)
round(prop.table(tabtp5miss), digits=3)

tfr$tp6miss <- as.numeric(is.na(tfr$TP6))
tabtp6miss <- table(tfr$tp6miss)
round(prop.table(tabtp6miss), digits=3)


tpmiss <- subset(tfr, select = c("tp1miss", "tp2miss", "tp3miss", "tp4miss", "tp5miss", "tp6miss"))
tfr$tpmiss <- rowSums(tpmiss)

tfr$tpmissout <- as.numeric(tfr$tpmiss==6)
tfr$tpmissout
tabtpmissout <- table(tfr$tpmissout)
prop.table(tabtpmissout)

##for OP
tfr$op1miss <- as.numeric(is.na(tfr$OP1))
tabop1miss <- table(tfr$op1miss)
round(prop.table(tabop1miss), digits=3)

tfr$op2miss <- as.numeric(is.na(tfr$OP2))
tabop2miss <- table(tfr$op2miss)
round(prop.table(tabop2miss), digits=3)

tfr$op3miss <- as.numeric(is.na(tfr$OP3))
tabop3miss <- table(tfr$op3miss)
round(prop.table(tabop3miss), digits=3)

tfr$op4miss <- as.numeric(is.na(tfr$OP4))
tabop4miss <- table(tfr$op4miss)
round(prop.table(tabop4miss), digits=3)

tfr$op5miss <- as.numeric(is.na(tfr$OP5))
tabop5miss <- table(tfr$op5miss)
round(prop.table(tabop5miss), digits=3)

tfr$op6miss <- as.numeric(is.na(tfr$OP6))
tabop6miss <- table(tfr$op6miss)
round(prop.table(tabop6miss), digits=3)

opmiss <- subset(tfr, select = c("op1miss", "op2miss", "op3miss", "op4miss", "op5miss", "op6miss"))
tfr$opmiss <- rowSums(opmiss)

tfr$opmissout <- as.numeric(tfr$opmiss==6)
tfr$opmissout
tabopmissout <- table(tfr$opmissout)
prop.table(tabopmissout)

##for SM

tfr$sm1miss <- as.numeric(is.na(tfr$SM1))
tabsm1miss <- table(tfr$sm1miss)
round(prop.table(tabsm1miss), digits=3)

tfr$sm2miss <- as.numeric(is.na(tfr$SM2))
tabsm2miss <- table(tfr$sm2miss)
round(prop.table(tabsm2miss), digits=3)

tfr$sm3miss <- as.numeric(is.na(tfr$SM3))
tabsm3miss <- table(tfr$sm3miss)
round(prop.table(tabsm3miss), digits=3)

tfr$sm4miss <- as.numeric(is.na(tfr$SM4))
tabsm4miss <- table(tfr$sm4miss)
round(prop.table(tabsm4miss), digits=3)
colnames(tfr)


smmiss <- subset(tfr, select = c("sm1miss", "sm2miss", "sm3miss", "sm4miss"))
tfr$smmiss <- rowSums(smmiss)

tfr$smmissout <- as.numeric(tfr$smmiss==12)
tfr$smmissout
tabsmmissout <- table(tfr$smmissout)
prop.table(tabsmmissout)

##for CM

tfr$cm1miss <- as.numeric(is.na(tfr$CM1))
tabcm1miss <- table(tfr$cm1miss)
round(prop.table(tabcm1miss), digits=3)


tfr$cm2miss <- as.numeric(is.na(tfr$CM2))
tabcm2miss <- table(tfr$cm2miss)
round(prop.table(tabcm2miss), digits=3)

tfr$cm3miss <- as.numeric(is.na(tfr$CM3))
tabcm3miss <- table(tfr$cm3miss)
round(prop.table(tabcm3miss), digits=3)

tfr$cm4miss <- as.numeric(is.na(tfr$CM4))
tabcm4miss <- table(tfr$cm4miss)
round(prop.table(tabcm4miss), digits=3)


cmmiss <- subset(tfr, select = c("cm1miss", "cm2miss", "cm3miss", "cm4miss"))
tfr$cmmiss <- rowSums(cmmiss)

tfr$cmmissout <- as.numeric(tfr$cmmiss ==4)
tfr$cmmissout
tabcmmissout <- table(tfr$cmmissout)
prop.table(tabcmmissout)

# calculate dropout rate due to missings

dropout <- subset(tfr, select = c("wmissout", "spmissout", "tpmissout", "opmissout", "smmissout", "cmmissout"))
tfr$dropout <- as.numeric(rowSums(dropout) == 6)
tabdropout <- table(tfr$dropout)
prop.table(tabdropout)

#delete cases with missings on all indicators for 1 causal factor

tfr = tfr[tfr$wmissout == 0,]
tfr = tfr[tfr$spmissout == 0,]
tfr = tfr[tfr$tpmissout == 0,]
tfr = tfr[tfr$opmissout == 0,]
tfr = tfr[tfr$smmissout == 0,]
tfr = tfr[tfr$cmmissout == 0,]
tfr

describe(tfr)

#check missings
missings <- as.numeric(complete.cases(tfr))
tabmiss <- table(missings)
prop.table(tabmiss)

#calculate dropout rate - from 1096 to 1087

do <- round(1-(1087/1096), digits=3)
do


#####CALIBRATION########

#calibrate indicator sets, recode neutral answers on the Likert scale into full nonmembership; calibrate indicator sets
#the values on the likert scales for powerlessness and meaninglessness need to be reverted because they should mean 
#powerfulness and meaningfulness
##direct method


tfr$wa <- tfr$W1
tfr$wa[tfr$W1 == 4] <- 3
tfr$wa[tfr$W1 == 5] <- 4
tfr$wa


tfr$wb <- tfr$W2
tfr$wb[tfr$W2 == 4] <- 3
tfr$wb[tfr$W2 == 5] <- 4

tfr$wc <- tfr$W3
tfr$wc[tfr$W3 == 4] <- 3
tfr$wc[tfr$W3 == 5] <- 4


tfr$wd <- tfr$W4
tfr$wd[tfr$W4 == 4] <- 3
tfr$wd[tfr$W4 == 5] <- 4


tfr$we <- tfr$W5
tfr$we[tfr$W5 == 4] <- 3
tfr$we[tfr$W5 == 5] <- 4


tfr$w1 <- calibrate(tfr$wa, type = "fuzzy", thresholds = c(1, 2.5, 4), logistic = TRUE)
tfr$w1[tfr$W1 == 3] <- 0.05
tfr$w1
tfr$w2 <- calibrate(tfr$wb, type = "fuzzy", thresholds = c(1, 2.5, 4), logistic = TRUE)
tfr$w2[tfr$W2 == 3] <- 0.05
tfr$w2
tfr$w3 <- calibrate(tfr$wc, type = "fuzzy", thresholds = c(1, 2.5, 4), logistic = TRUE)
tfr$w3[tfr$W3 == 3] <- 0.05
tfr$w3
tfr$w4 <- calibrate(tfr$wd, type = "fuzzy", thresholds = c(1, 2.5, 4), logistic = TRUE)
tfr$w4[tfr$W4 == 3] <- 0.05
tfr$w4
tfr$w5 <- calibrate(tfr$we, type = "fuzzy", thresholds = c(1, 2.5, 4), logistic = TRUE)
tfr$w5[tfr$W5 == 3] <- 0.05
tfr$w5

## for SP
tfr$spa <- tfr$SP1
tfr$spa[tfr$SP1 == 4] <- 3
tfr$spa[tfr$SP1 == 5] <- 4
tfr$spa


tfr$spb <- tfr$SP2
tfr$spb[tfr$SP2 == 4] <- 3
tfr$spb[tfr$SP2 == 5] <- 4

tfr$spc <- tfr$SP3
tfr$spc[tfr$SP3 == 4] <- 3
tfr$spc[tfr$SP3 == 5] <- 4


tfr$spd <- tfr$SP4
tfr$spd[tfr$SP4 == 4] <- 3
tfr$spd[tfr$SP4 == 5] <- 4


tfr$spe <- tfr$SP5
tfr$spe[tfr$SP5 == 4] <- 3
tfr$spe[tfr$SP5 == 5] <- 4

tfr$spf <- tfr$SP6
tfr$spf[tfr$SP6 == 4] <- 3
tfr$spf[tfr$SP6 == 5] <- 4


tfr$sp1 <- calibrate(tfr$spa, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sp1[tfr$SP1 == 3] <- 0.05
tfr$sp1
tfr$sp2 <- calibrate(tfr$spb, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sp2[tfr$SP2 == 3] <- 0.05
tfr$sp2
tfr$sp3 <- calibrate(tfr$spc, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sp3[tfr$SP3 == 3] <- 0.05
tfr$sp3
tfr$sp4 <- calibrate(tfr$spd, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sp4[tfr$SP4 == 3] <- 0.05
tfr$sp4
tfr$sp5 <- calibrate(tfr$spe, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sp5[tfr$SP5 == 3] <- 0.05
tfr$sp5
tfr$sp6 <- calibrate(tfr$spf, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sp6[tfr$SP6 == 3] <- 0.05
tfr$sp6

## for TP
tfr$tpa <- tfr$TP1
tfr$tpa[tfr$TP1 == 4] <- 3
tfr$tpa[tfr$TP1 == 5] <- 4
tfr$tpa


tfr$tpb <- tfr$TP2
tfr$tpb[tfr$TP2 == 4] <- 3
tfr$tpb[tfr$TP2 == 5] <- 4

tfr$tpc <- tfr$TP3
tfr$tpc[tfr$TP3 == 4] <- 3
tfr$tpc[tfr$TP3 == 5] <- 4


tfr$tpd <- tfr$TP4
tfr$tpd[tfr$TP4 == 4] <- 3
tfr$tpd[tfr$TP4 == 5] <- 4


tfr$tpe <- tfr$TP5
tfr$tpe[tfr$TP5 == 4] <- 3
tfr$tpe[tfr$TP5 == 5] <- 4

tfr$tpf <- tfr$TP6
tfr$tpf[tfr$TP6 == 4] <- 3
tfr$tpf[tfr$TP6 == 5] <- 4


tfr$tp1 <- calibrate(tfr$tpa, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$tp1[tfr$TP1 == 3] <- 0.05
tfr$tp1
tfr$tp2 <- calibrate(tfr$tpb, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$tp2[tfr$TP2 == 3] <- 0.05
tfr$tp2
tfr$tp3 <- calibrate(tfr$tpc, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$tp3[tfr$TP3 == 3] <- 0.05
tfr$tp3
tfr$tp4 <- calibrate(tfr$tpd, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$tp4[tfr$TP4 == 3] <- 0.05
tfr$tp4
tfr$tp5 <- calibrate(tfr$tpe, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$tp5[tfr$TP5 == 3] <- 0.05
tfr$tp5
tfr$tp6 <- calibrate(tfr$tpf, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$tp6[tfr$TP6 == 3] <- 0.05
tfr$tp6

##for OP

tfr$opa <- tfr$OP1
tfr$opa[tfr$OP1 == 4] <- 3
tfr$opa[tfr$OP1 == 5] <- 4
tfr$opa


tfr$opb <- tfr$OP2
tfr$opb[tfr$OP2 == 4] <- 3
tfr$opb[tfr$OP2 == 5] <- 4

tfr$opc <- tfr$OP3
tfr$opc[tfr$OP3 == 4] <- 3
tfr$opc[tfr$OP3 == 5] <- 4


tfr$opd <- tfr$OP4
tfr$opd[tfr$OP4 == 4] <- 3
tfr$opd[tfr$OP4 == 5] <- 4


tfr$ope <- tfr$OP5
tfr$ope[tfr$OP5 == 4] <- 3
tfr$ope[tfr$OP5 == 5] <- 4

tfr$opf <- tfr$OP6
tfr$opf[tfr$OP6 == 4] <- 3
tfr$opf[tfr$OP6 == 5] <- 4


tfr$op1 <- calibrate(tfr$opa, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$op1[tfr$OP1 == 3] <- 0.05
tfr$op1
tfr$op2 <- calibrate(tfr$opb, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$op2[tfr$OP2 == 3] <- 0.05
tfr$op2
tfr$op3 <- calibrate(tfr$opc, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$op3[tfr$OP3 == 3] <- 0.05
tfr$op3
tfr$op4 <- calibrate(tfr$opd, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$op4[tfr$OP4 == 3] <- 0.05
tfr$op4
tfr$op5 <- calibrate(tfr$ope, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$op5[tfr$OP5 == 3] <- 0.05
tfr$op5
tfr$op6 <- calibrate(tfr$opf, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$op6[tfr$OP6 == 3] <- 0.05
tfr$op6

##for SM

tfr$sma <- tfr$SM1
tfr$sma[tfr$SM1 == 4] <- 3
tfr$sma[tfr$SM1 == 5] <- 4
tfr$sma


tfr$smb <- tfr$SM2
tfr$smb[tfr$SM2 == 4] <- 3
tfr$smb[tfr$SM2 == 5] <- 4

tfr$smc <- tfr$SM3
tfr$smc[tfr$SM3 == 4] <- 3
tfr$smc[tfr$SM3 == 5] <- 4


tfr$smd <- tfr$SM4
tfr$smd[tfr$SM4 == 4] <- 3
tfr$smd[tfr$SM4 == 5] <- 4

colnames(tfr)



tfr$sm1 <- calibrate(tfr$sma, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sm1[tfr$SM1 == 3] <- 0.05
tfr$sm1
tfr$sm2 <- calibrate(tfr$smb, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sm2[tfr$SM2 == 3] <- 0.05
tfr$sm2
tfr$sm3 <- calibrate(tfr$smc, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sm3[tfr$SM3 == 3] <- 0.05
tfr$sm3
tfr$sm4 <- calibrate(tfr$smd, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sm4[tfr$SM4 == 3] <- 0.05
tfr$sm4


##For CM
## for CM
tfr$cma <- tfr$CM1
tfr$cma[tfr$CM1 == 4] <- 3
tfr$cma[tfr$CM1 == 5] <- 4
tfr$cma


tfr$cmb <- tfr$CM2
tfr$cmb[tfr$CM2 == 4] <- 3
tfr$cmb[tfr$CM2 == 5] <- 4

tfr$cmc <- tfr$CM3
tfr$cmc[tfr$CM3 == 4] <- 3
tfr$cmc[tfr$CM3 == 5] <- 4


tfr$cmd <- tfr$CM4
tfr$cmd[tfr$CM4 == 4] <- 3
tfr$cmd[tfr$CM4 == 5] <- 4

tfr$cm1 <- calibrate(tfr$cma, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$cm1[tfr$CM1 == 3] <- 0.05
tfr$cm1
tfr$cm2 <- calibrate(tfr$cmb, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$cm2[tfr$CM2 == 3] <- 0.05
tfr$cm2
tfr$cm3 <- calibrate(tfr$cmc, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$cm3[tfr$CM3 == 3] <- 0.05
tfr$cm3
tfr$cm4 <- calibrate(tfr$cmd, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$cm4[tfr$CM4 == 3] <- 0.05
tfr$cm4

skewcm1 <- as.numeric(tfr$cm1 > 0.5)
tabcm1 <- table(skewcm1)
prop.table(tabcm1)

skewcm2 <- as.numeric(tfr$cm2 > 0.5)
tabcm2 <- table(skewcm2)
prop.table(tabcm2)

skewcm3 <- as.numeric(tfr$cm3 > 0.5)
tabcm3 <- table(skewcm3)
prop.table(tabcm3)

skewcm4 <- as.numeric(tfr$cm4 > 0.5)
tabcm4 <- table(skewcm4)
prop.table(tabcm4)

#####Aggregation########

#build aggregated sets using the maximum rule, check for values of 0.5 and for skewedness

tfr$w <- pmax(tfr$w1, tfr$w2, tfr$w3, tfr$w4, tfr$w5, na.rm = TRUE)
tfr$w

checkw <- as.numeric(tfr$w == 0.5)
sum(checkw == 1)

wskew <- as.numeric(tfr$w > 0.5)
tabw <- table(wskew)
prop.table(tabw)

tfr$sp <- pmax(tfr$sp1, tfr$sp2, tfr$sp3, tfr$sp4, tfr$sp5, tfr$sp6, na.rm = TRUE)
tfr$sp

checksp <- as.numeric(tfr$sp == 0.5)
sum(checksp == 1)

skewsp <- as.numeric(tfr$sp > 0.5)
tabsp <- table(skewsp)
prop.table(tabsp)


tfr$tp <- pmax(tfr$tp1, tfr$tp2, tfr$tp3, tfr$tp4, tfr$tp5, tfr$tp6, na.rm = TRUE)
tfr$tp

checktp <- as.numeric(tfr$tp == 0.5)
sum(checktp == 1)

skewtp <- as.numeric(tfr$tp > 0.5)
tabtp <- table(skewtp)
prop.table(tabtp)

tfr$op <- pmax(tfr$op1, tfr$op2, tfr$op3, tfr$op4, tfr$op5, tfr$op6, na.rm = TRUE)
tfr$op

checkop <- as.numeric(tfr$op == 0.5)
sum(checkop == 1)

skewop <- as.numeric(tfr$op > 0.5)
tabop <- table(skewop)
prop.table(tabop)

tfr$sm <- pmax(tfr$sm1, tfr$sm2, tfr$sm3, tfr$sm4, na.rm = TRUE)
tfr$sm

checksm <- as.numeric(tfr$sm == 0.5)
sum(checksm == 1)

skewsm <- as.numeric(tfr$sm > 0.5)
tabsm <- table(skewsm)
prop.table(tabsm)


tfr$cm <- pmax(tfr$cm1, tfr$cm2, tfr$cm3, tfr$cm4, na.rm = TRUE)
tfr$cm

checkcm <- as.numeric(tfr$cm == 0.5)
sum(checkcm == 1)

skewcm <- as.numeric(tfr$cm > 0.5)
tabcm <- table(skewcm)
prop.table(tabcm)

#exclude cases with set membership 0.5

tfr = tfr[tfr$w != 0.5,]
tfr = tfr[tfr$op != 0.5,] 
tfr = tfr[tfr$sp != 0.5,]
tfr = tfr[tfr$tp != 0.5,]
tfr = tfr[tfr$sm != 0.5,]
tfr = tfr[tfr$cm != 0.5,]
tfr
describe(tfr)
colnames(tfr)

# rename raw variables
tfr <- rename(tfr, W = w)
tfr <- rename(tfr, SP=  sp)
tfr <- rename(tfr, TP = tp)
tfr <- rename(tfr, OP = op)
tfr <- rename(tfr, SM = sm)
tfr <- rename(tfr, CM = cm)
colnames (tfr)


tff_d2_c1 <- subset(tfr, select = c("W", "SP", "TP", "OP", "SM", "CM"))
write.csv2(tff_d2_c1, file = "tff_d2_c1.csv")



###################ANALYSIS OF NECESSITY############

tff <- read.csv("tff_d2_c1.csv", row.names=1, header = TRUE, sep = ";",  dec = ",")
describe(tff)

#analysis of necessity for W and w  

superSubset(tff, outcome = "W", 
            conditions = "OP, SP, TP, SM, CM", 
            incl.cut = 0.9, cov.cut = 0.6)

superSubset(tff, outcome = "~W", 
            conditions = "OP, SP, TP, SM, CM", 
            incl.cut = 0.9, cov.cut = 0.6)

ttw1


#calculate Z-Test for necessary conditions. I apply a RoN threshold of 0.55, so only for NCs that meet this threshold

nw <- 1-tff$W

# for NC1 TP+cm

NC1W <- compute("TP+cm", data=tff)
obsNC1W <- as.numeric(NC1W >= tff$W)

zNC1W <- round(((((sum(obsNC1W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC1W

pzNC1W <- 2*pnorm(-abs(zNC1W))
pzNC1W

# for NC2 TP+sm

NC2W <- compute("TP+sm", data=tff)
obsNC2W <- as.numeric(NC2W >= tff$W)

zNC2W <- round(((((sum(obsNC2W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC2W

pzNC2W <- 2*pnorm(-abs(zNC2W))
pzNC2W

# for NC3 OP+cm

NC3W <- compute("OP+cm", data=tff)
obsNC3W <- as.numeric(NC3W >= tff$W)

zNC3W <- round(((((sum(obsNC3W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC3W

pzNC3W <- 2*pnorm(-abs(zNC3W))
pzNC3W

# for NC4 tp+SM+cm

NC4W <- compute("tp+SM+cm", data=tff)
obsNC4W <- as.numeric(NC4W >= tff$W)

zNC4W <- round(((((sum(obsNC4W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC4W

pzNC4W <- 2*pnorm(-abs(zNC4W))
pzNC4W

# for NC5 sp+TP+SM

NC5W <- compute("sp+TP+SM", data=tff)
obsNC5W <- as.numeric(NC5W >= tff$W)

zNC5W <- round(((((sum(obsNC5W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC5W

pzNC5W <- 2*pnorm(-abs(zNC5W))
pzNC5W

# for NC6 SP+sm+cm

NC6W <- compute("SP+sm+cm", data=tff)
obsNC6W <- as.numeric(NC6W >= tff$W)

zNC6W <- round(((((sum(obsNC6W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC6W

pzNC6W <- 2*pnorm(-abs(zNC6W))
pzNC6W

# for NC7 SP+sm+CM

NC7W <- compute("SP+sm+CM", data=tff)
obsNC7W <- as.numeric(NC7W >= tff$W)

zNC7W <- round(((((sum(obsNC7W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC7W

pzNC7W <- 2*pnorm(-abs(zNC7W))
pzNC7W

# for NC9 SP+tp+cm

NC9W <- compute("SP+tp+cm", data=tff)
obsNC9W <- as.numeric(NC9W >= tff$W)

zNC9W <- round(((((sum(obsNC9W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC9W

pzNC9W <- 2*pnorm(-abs(zNC9W))
pzNC9W

# for NC10 SP+TP+CM

NC10W <- compute("SP+TP+CM", data=tff)
obsNC10W <- as.numeric(NC10W >= tff$W)

zNC10W <- round(((((sum(obsNC10W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC10W

pzNC10W <- 2*pnorm(-abs(zNC10W))
pzNC10W

# for NC11 SP+TP+SM

NC11W <- compute("SP+TP+SM", data=tff)
obsNC11W <- as.numeric(NC11W >= tff$W)

zNC11W <- round(((((sum(obsNC11W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC11W

pzNC11W <- 2*pnorm(-abs(zNC11W))
pzNC11W

# for NC12 op+TP+CM

NC12W <- compute("op+TP+CM", data=tff)
obsNC12W <- as.numeric(NC12W >= tff$W)

zNC12W <- round(((((sum(obsNC12W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC12W

pzNC12W <- 2*pnorm(-abs(zNC12W))
pzNC12W

# for NC13 op+TP+SM

NC13W <- compute("op+TP+SM", data=tff)
obsNC13W <- as.numeric(NC13W >= tff$W)

zNC13W <- round(((((sum(obsNC13W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC13W

pzNC13W <- 2*pnorm(-abs(zNC13W))
pzNC13W

# for NC14 op+sp+TP

NC14W <- compute("op+sp+TP", data=tff)
obsNC14W <- as.numeric(NC14W >= tff$W)

zNC14W <- round(((((sum(obsNC14W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC14W

pzNC14W <- 2*pnorm(-abs(zNC14W))
pzNC14W

# for NC15 op+SP+cm

NC15W <- compute("op+SP+cm", data=tff)
obsNC15W <- as.numeric(NC15W >= tff$W)

zNC15W <- round(((((sum(obsNC15W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC15W

pzNC15W <- 2*pnorm(-abs(zNC15W))
pzNC15W

# for NC16 op+SP+TP

NC16W <- compute("op+SP+TP", data=tff)
obsNC16W <- as.numeric(NC16W >= tff$W)

zNC16W <- round(((((sum(obsNC16W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC16W

pzNC16W <- 2*pnorm(-abs(zNC16W))
pzNC16W

# for NC17 OP+sm+CM

NC17W <- compute("OP+sm+CM", data=tff)
obsNC17W <- as.numeric(NC17W >= tff$W)

zNC17W <- round(((((sum(obsNC17W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC17W

pzNC17W <- 2*pnorm(-abs(zNC17W))
pzNC17W

# for NC18 OP+tp+sm

NC18W <- compute("OP+tp+sm", data=tff)
obsNC18W <- as.numeric(NC18W >= tff$W)

zNC18W <- round(((((sum(obsNC18W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC18W

pzNC18W <- 2*pnorm(-abs(zNC18W))
pzNC18W

# for NC19 OP+TP+SM

NC19W <- compute("OP+TP+SM", data=tff)
obsNC19W <- as.numeric(NC19W >= tff$W)

zNC19W <- round(((((sum(obsNC19W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC19W

pzNC19W <- 2*pnorm(-abs(zNC19W))
pzNC19W

# for NC20 OP+sp+sm

NC20W <- compute("OP+sp+sm", data=tff)
obsNC20W <- as.numeric(NC20W >= tff$W)

zNC20W <- round(((((sum(obsNC20W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC20W

pzNC20W <- 2*pnorm(-abs(zNC20W))
pzNC20W

# for NC21 OP+sp+TP

NC21W <- compute("OP+sp+TP", data=tff)
obsNC21W <- as.numeric(NC21W >= tff$W)

zNC21W <- round(((((sum(obsNC21W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC21W

pzNC21W <- 2*pnorm(-abs(zNC21W))
pzNC21W

# for NC22 OP+SP+SM

NC22W <- compute("OP+SP+SM", data=tff)
obsNC22W <- as.numeric(NC22W >= tff$W)

zNC22W <- round(((((sum(obsNC22W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC22W

pzNC22W <- 2*pnorm(-abs(zNC22W))
pzNC22W

# for NC24 OP+SP+tp

NC24W <- compute("OP+SP+tp", data=tff)
obsNC24W <- as.numeric(NC24W >= tff$W)

zNC24W <- round(((((sum(obsNC24W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC24W

pzNC24W <- 2*pnorm(-abs(zNC24W))
pzNC24W

# for NC25 OP+SP+TP

NC25W <- compute("OP+SP+TP", data=tff)
obsNC25W <- as.numeric(NC25W >= tff$W)

zNC25W <- round(((((sum(obsNC25W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC25W

pzNC25W <- 2*pnorm(-abs(zNC25W))
pzNC25W

# for NC26 op+sp+SM+cm

NC26W <- compute("op+sp+SM+cm ", data=tff)
obsNC26W <- as.numeric(NC26W >= tff$W)

zNC26W <- round(((((sum(obsNC26W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC26W

pzNC26W <- 2*pnorm(-abs(zNC26W))
pzNC26W

# for NC27 op+SP+tp+sm

NC27W <- compute("op+SP+tp+sm", data=tff)
obsNC27W <- as.numeric(NC27W >= tff$W)

zNC27W <- round(((((sum(obsNC27W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC27W

pzNC27W <- 2*pnorm(-abs(zNC27W))
pzNC27W

# for NC28 OP+sp+tp+SM

NC28W <- compute("OP+sp+tp+SM", data=tff)
obsNC28W <- as.numeric(NC28W >= tff$W)

zNC28W <- round(((((sum(obsNC28W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC28W

pzNC28W <- 2*pnorm(-abs(zNC28W))
pzNC28W

# for NC29 op+SP+tp+SM+CM

NC29W <- compute("op+SP+tp+SM+CM", data=tff)
obsNC29W <- as.numeric(NC29W >= tff$W)

zNC29W <- round(((((sum(obsNC29W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC29W

pzNC29W <- 2*pnorm(-abs(zNC29W))
pzNC29W





#Test for skewness of necessary conditions. It is tested whether the proportion of cases in all but one 
#specific area of the XY plot is lower than ten per cent.
##First, visual check with XY Plots; then calculate number of cases outside each area; 
#and calulcate percentage outside area; 
#there are a lot of cases on the diagonal, which are attributed to two area. This is why 
#the percentages don't add up to 100.


##for SP + TP + OP

A1NC25W <- round(sum(as.numeric((NC25W <= tff$W) & (NC25W >= nw)))/(1087/100), digits = 1)

A2NC25W <- round(sum(as.numeric((NC25W >= tff$W)&(NC25W >= nw)))/(1087/100), digits = 1)

A3NC25W <- round(sum(as.numeric((NC25W >= tff$W)&(NC25W <= nw)))/(1087/100), digits=1)

A4NC25W <- round(sum(as.numeric((NC25W <= tff$W)&(NC25W <= nw)))/(1087/100), digits = 1)

sum(A1NC25W, A2NC25W, A3NC25W, A4NC25W)

par(mfrow=c(3, 2))
plot(tff$W, NC25W, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness test for SP + TP + OP ',
     xlab=' SP + TP + OP ',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1NC25W)
text(x=0.8, y=0.5, labels = A2NC25W)
text(x=0.5, y=0.2, labels = A3NC25W)
text(x=0.2, y=0.5, labels = A4NC25W)


#################ANALYSIS OF SUFFICIENCY FOR W##################

tff <- read.csv("tff_d2_c1.csv", row.names=1, header = TRUE, sep = ";",  dec = ",")
describe(tff)

nw <- 1-tff$W

#analysis of sufficiency for W
# calculate untenable assumptions

nec <- deMorgan("SP + TP + OP", snames = "SP, TP, OP")
nec

#### W1: frequency threshold 5, raw consistency threshold 0.937 ####

ttW1 <- truthTable(tff, outcome="W", 
                        conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.937, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                        complete=FALSE, show.cases=FALSE)
ttW1


#parsimonious solution
psW1 <- eqmcc(ttW1, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW1

#exclude contradictory assumptions & rows: sp*tp*op (contradict statement of necessity)

ttW1new <- ttW1
nec
rows <- findRows("sp*tp*op", ttW1new, remainders = FALSE)
rows

ttW1new$tt[as.character(rows), "OUT"] <- 0
ttW1new


#conservative solution
csW1 <- eqmcc(ttW1new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW1

#enhanced parsimonious solution (all simplifying assumptions)

epsW1 <- eqmcc(ttW1new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
              rowdom=TRUE, use.tilde=FALSE)
epsW1


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  SP*SM*CM + TP*OP*SM*CM

epsolW1 <- compute("SP*SM*CM + TP*OP*SM*CM", data=tff)


obsepsolW1 <- as.numeric(epsolW1 <= tff$W)
zepsolW1 <- round(((((sum(obsepsolW1 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zepsolW1
pzepsolW1 <- 2*pnorm(-abs(zepsolW1))
pzepsolW1


#intermediate solution with directional expectations SM -> W, CM -> W
isW1 <- eqmcc(ttW1new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW1


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: SP*SM*CM + TP*OP*SM*CM

path1W1 <- compute("SP*SM*CM", data=tff)
path2W1 <- compute("TP*OP*SM*CM", data=tff)


A1path1W1 <- round(sum(as.numeric((path1W1 <= tff$W) & (path1W1 >= nw)))/(1087/100), digits = 1)

A2path1W1 <- round(sum(as.numeric((path1W1 >= tff$W)&(path1W1 >= nw)))/(1087/100), digits = 1)

A3path1W1 <- round(sum(as.numeric((path1W1 >= tff$W)&(path1W1 <= nw)))/(1087/100), digits=1)

A4path1W1 <- round(sum(as.numeric((path1W1 <= tff$W)&(path1W1 <= nw)))/(1087/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W1)
text(x=0.8, y=0.5, labels = A2path1W1)
text(x=0.5, y=0.2, labels = A3path1W1)
text(x=0.2, y=0.5, labels = A4path1W1)


A1path2W1 <- round(sum(as.numeric((path2W1 <= tff$W) & (path2W1 >= nw)))/(1087/100), digits = 1)

A2path2W1 <- round(sum(as.numeric((path2W1 >= tff$W)&(path2W1 >= nw)))/(1087/100), digits = 1)

A3path2W1 <- round(sum(as.numeric((path2W1 >= tff$W)&(path2W1 <= nw)))/(1087/100), digits=1)

A4path2W1 <- round(sum(as.numeric((path2W1 <= tff$W)&(path2W1 <= nw)))/(1087/100), digits = 1)


plot(tff$W, path2W1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab='Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W1)
text(x=0.8, y=0.5, labels = A2path2W1)
text(x=0.5, y=0.2, labels = A3path2W1)
text(x=0.2, y=0.5, labels = A4path2W1)



A1epsolW1 <- round(sum(as.numeric((epsolW1 <= tff$W) & (epsolW1 >= nw)))/(1087/100), digits = 1)

A2epsolW1 <- round(sum(as.numeric((epsolW1 >= tff$W)&(epsolW1 >= nw)))/(1087/100), digits = 1)

A3epsolW1 <- round(sum(as.numeric((epsolW1 >= tff$W)&(epsolW1 <= nw)))/(1087/100), digits=1)

A4epsolW1 <- round(sum(as.numeric((epsolW1 <= tff$W)&(epsolW1 <= nw)))/(1087/100), digits = 1)


isW1
plot(tff$W, epsolW1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW1)
text(x=0.8, y=0.5, labels = A2epsolW1)
text(x=0.5, y=0.2, labels = A3epsolW1)
text(x=0.2, y=0.5, labels = A4epsolW1)

#####W2  frequency threshold 5, raw consistency threshold 0.934 #####


ttW2 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.934, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW2

#parsimonious solution
psW2 <- eqmcc(ttW2, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW2

#exclude contradictory assumptions & rows: sp*tp*op (contradict statement of necessity)

ttW2new <- ttW2
nec
rows <- findRows("sp*tp*op", ttW2new, remainders = FALSE)
rows

ttW2new$tt[as.character(rows), "OUT"] <- 0
ttW2new


#conservative solution
csW2 <- eqmcc(ttW2new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW2

#enhanced parsimonious solution (all simplifying assumptions)

epsW2 <- eqmcc(ttW2new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW2


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  OP*SM*CM + SP*SM*CM + TP*SM*CM
epsolW2 <- compute("OP*SM*CM + SP*SM*CM + TP*SM*CM ", data=tff)

obsepsolW2 <- as.numeric(epsolW2 <= tff$W)

zepsolW2 <- round(((((sum(obsepsolW2 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zepsolW2
pzepsolW2 <- 2*pnorm(-abs(zepsolW2))
pzepsolW2


#intermediate solution with directional expectations SM -> W, CM -> W
isW2 <- eqmcc(ttW2new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW2

epsW2

# identify simplifying assumptions and easy counterfactuals

SAW2 <- psW2$SA
SAW2

ECW2 <- isW2$i.sol$C1P1$EC
ECW2


#Z-Test for intermediate solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  OP*SM*CM + SP*SM*CM + TP*SM*CM
isW2
isolW2 <- compute("OP*SM*CM + SP*SM*CM + TP*SM*CM ", data=tff)

obsisolW2 <- as.numeric(isolW2 <= tff$W)
zisolW2 <- round(((((sum(obsisolW2 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zisolW2
pzisolW2 <- 2*pnorm(-abs(zisolW2))
pzisolW2

# for path 1
# Calculate membership in path  OP*SM*CM
path1isolW2 <- compute("OP*SM*CM ", data=tff)

obspath1isolW2 <- as.numeric(path1isolW2 <= tff$W)
zpath1isolW2 <- round(((((sum(obspath1isolW2 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zpath1isolW2
pzpath1isolW2 <- 2*pnorm(-abs(zpath1isolW2))
pzpath1isolW2

#for path 2
# Calculate membership in path  SP*SM*CM
path2isolW2 <- compute("SP*SM*CM ", data=tff)

obspath2isolW2 <- as.numeric(path2isolW2 <= tff$W)
zpath2isolW2 <- round(((((sum(obspath2isolW2 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zpath2isolW2
pzpath2isolW2 <- 2*pnorm(-abs(zpath2isolW2))
pzpath2isolW2

# for path 3
# Calculate membership in path  TP*SM*CM
path3isolW2 <- compute("TP*SM*CM ", data=tff)

obspath3isolW2 <- as.numeric(path3isolW2 <= tff$W)
zpath3isolW2 <- round(((((sum(obspath3isolW2 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zpath3isolW2
pzpath3isolW2 <- 2*pnorm(-abs(zpath3isolW2))
pzpath3isolW2



#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: OP*SM*CM + SP*SM*CM + TP*SM*CM

path1W2 <- compute("OP*SM*CM ", data=tff)
path2W2 <- compute("SP*SM*CM ", data=tff)
path3W2 <- compute("TP*SM*CM ", data=tff)


A1path1W2 <- round(sum(as.numeric((path1W2 <= tff$W) & (path1W2 >= nw)))/(1087/100), digits = 1)

A2path1W2 <- round(sum(as.numeric((path1W2 >= tff$W)&(path1W2 >= nw)))/(1087/100), digits = 1)

A3path1W2 <- round(sum(as.numeric((path1W2 >= tff$W)&(path1W2 <= nw)))/(1087/100), digits=1)

A4path1W2 <- round(sum(as.numeric((path1W2 <= tff$W)&(path1W2 <= nw)))/(1087/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W2)
text(x=0.8, y=0.5, labels = A2path1W2)
text(x=0.5, y=0.2, labels = A3path1W2)
text(x=0.2, y=0.5, labels = A4path1W2)



A1path2W2 <- round(sum(as.numeric((path2W2 <= tff$W) & (path2W2 >= nw)))/(1087/100), digits = 1)

A2path2W2 <- round(sum(as.numeric((path2W2 >= tff$W)&(path2W2 >= nw)))/(1087/100), digits = 1)

A3path2W2 <- round(sum(as.numeric((path2W2 >= tff$W)&(path2W2 <= nw)))/(1087/100), digits=1)

A4path2W2 <- round(sum(as.numeric((path2W2 <= tff$W)&(path2W2 <= nw)))/(1087/100), digits = 1)

plot(tff$W, path2W2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab='Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W2)
text(x=0.8, y=0.5, labels = A2path2W2)
text(x=0.5, y=0.2, labels = A3path2W2)
text(x=0.2, y=0.5, labels = A4path2W2)

A1path3W2 <- round(sum(as.numeric((path3W2 <= tff$W) & (path3W2 >= nw)))/(1087/100), digits = 1)

A2path3W2 <- round(sum(as.numeric((path3W2 >= tff$W)&(path3W2 >= nw)))/(1087/100), digits = 1)

A3path3W2 <- round(sum(as.numeric((path3W2 >= tff$W)&(path3W2 <= nw)))/(1087/100), digits=1)

A4path3W2 <- round(sum(as.numeric((path3W2 <= tff$W)&(path3W2 <= nw)))/(1087/100), digits = 1)

plot(tff$W, path3W2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness path 3',
     xlab='Path 3',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3W2)
text(x=0.8, y=0.5, labels = A2path3W2)
text(x=0.5, y=0.2, labels = A3path3W2)
text(x=0.2, y=0.5, labels = A4path3W2)


A1epsolW2 <- round(sum(as.numeric((epsolW2 <= tff$W) & (epsolW2 >= nw)))/(1087/100), digits = 1)

A2epsolW2 <- round(sum(as.numeric((epsolW2 >= tff$W)&(epsolW2 >= nw)))/(1087/100), digits = 1)

A3epsolW2 <- round(sum(as.numeric((epsolW2 >= tff$W)&(epsolW2 <= nw)))/(1087/100), digits=1)

A4epsolW2 <- round(sum(as.numeric((epsolW2 <= tff$W)&(epsolW2 <= nw)))/(1087/100), digits = 1)

isW2
plot(tff$W, epsolW2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW2)
text(x=0.8, y=0.5, labels = A2epsolW2)
text(x=0.5, y=0.2, labels = A3epsolW2)
text(x=0.2, y=0.5, labels = A4epsolW2)



####W3: frequency threshold 5, raw consistency threshold 0.932#######

ttW3 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.932, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW3


#parsimonious solution
psW3 <- eqmcc(ttW3, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW3

#exclude contradictory assumptions & rows: sptpop (contradict statement of necessity)

ttW3new <- ttW3
nec
rows <- findRows("sp*tp*op", ttW3new, remainders = FALSE)
rows

ttW3new$tt[as.character(rows), "OUT"] <- 0
ttW3new


#conservative solution
csW3 <- eqmcc(ttW3new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW3

#enhanced parsimonious solution (all simplifying assumptions)

epsW3 <- eqmcc(ttW3new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW3

#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term OP*SM*CM + TP*SM*CM + SP*tp*op*CM
epsolW3 <- compute("OP*SM*CM + TP*SM*CM + SP*tp*op*CM", data=tff)

obsepsolW3 <- as.numeric(epsolW3 <= tff$W)
zepsolW3 <- round(((((sum(obsepsolW3 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zepsolW3
pzepsolW3 <- 2*pnorm(-abs(zepsolW3))
pzepsolW3



#intermediate solution with directional expectations SM -> W, CM -> W
isW3 <- eqmcc(ttW3new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW3

# identify simplifying assumptions and easy counterfactuals

SAW3 <- psW3$SA
SAW3

ECW3 <- isW3$i.sol$C1P1$EC
ECW3


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: OP*SM*CM + TP*SM*CM + SP*tp*op*CM

path1W3 <- compute("OP*SM*CM", data=tff)
path2W3 <- compute("TP*SM*CM ", data=tff)
path3W3 <- compute("SP*tp*op*CM ", data=tff)

A1path1W3 <- round(sum(as.numeric((path1W3 <= tff$W) & (path1W3 >= nw)))/(1087/100), digits = 1)

A2path1W3 <- round(sum(as.numeric((path1W3 >= tff$W)&(path1W3 >= nw)))/(1087/100), digits = 1)

A3path1W3 <- round(sum(as.numeric((path1W3 >= tff$W)&(path1W3 <= nw)))/(1087/100), digits=1)

A4path1W3 <- round(sum(as.numeric((path1W3 <= tff$W)&(path1W3 <= nw)))/(1087/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W3)
text(x=0.8, y=0.5, labels = A2path1W3)
text(x=0.5, y=0.2, labels = A3path1W3)
text(x=0.2, y=0.5, labels = A4path1W3)



A1path2W3 <- round(sum(as.numeric((path2W3 <= tff$W) & (path2W3 >= nw)))/(1087/100), digits = 1)

A2path2W3 <- round(sum(as.numeric((path2W3 >= tff$W)&(path2W3 >= nw)))/(1087/100), digits = 1)

A3path2W3 <- round(sum(as.numeric((path2W3 >= tff$W)&(path2W3 <= nw)))/(1087/100), digits=1)

A4path2W3 <- round(sum(as.numeric((path2W3 <= tff$W)&(path2W3 <= nw)))/(1087/100), digits = 1)

plot(tff$W, path2W3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab='Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W3)
text(x=0.8, y=0.5, labels = A2path2W3)
text(x=0.5, y=0.2, labels = A3path2W3)
text(x=0.2, y=0.5, labels = A4path2W3)

A1path3W3 <- round(sum(as.numeric((path3W3 <= tff$W) & (path3W3 >= nw)))/(1087/100), digits = 1)

A2path3W3 <- round(sum(as.numeric((path3W3 >= tff$W)&(path3W3 >= nw)))/(1087/100), digits = 1)

A3path3W3 <- round(sum(as.numeric((path3W3 >= tff$W)&(path3W3 <= nw)))/(1087/100), digits=1)

A4path3W3 <- round(sum(as.numeric((path3W3 <= tff$W)&(path3W3 <= nw)))/(1087/100), digits = 1)

plot(tff$W, path3W3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness path 3',
     xlab='Path 3',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3W3)
text(x=0.8, y=0.5, labels = A2path3W3)
text(x=0.5, y=0.2, labels = A3path3W3)
text(x=0.2, y=0.5, labels = A4path3W3)

A1epsolW3 <- round(sum(as.numeric((epsolW3 <= tff$W) & (epsolW3 >= nw)))/(1087/100), digits = 1)

A2epsolW3 <- round(sum(as.numeric((epsolW3 >= tff$W)&(epsolW3 >= nw)))/(1087/100), digits = 1)

A3epsolW3 <- round(sum(as.numeric((epsolW3 >= tff$W)&(epsolW3 <= nw)))/(1087/100), digits=1)

A4epsolW3 <- round(sum(as.numeric((epsolW3 <= tff$W)&(epsolW3 <= nw)))/(1087/100), digits = 1)

isW3
plot(tff$W, epsolW3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW3)
text(x=0.8, y=0.5, labels = A2epsolW3)
text(x=0.5, y=0.2, labels = A3epsolW3)
text(x=0.2, y=0.5, labels = A4epsolW3)


#Z-Test for intermediate solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  OP*CM + SP*TP*CM + sp*TP*OP*SM + SP*TP*op*SM
isW3
isolW3 <- compute("OP*SM*CM + TP*SM*CM + SP*tp*op*CM", data=tff)

obsisolW3 <- as.numeric(isolW3 <= tff$W)
zisolW3 <- round(((((sum(obsisolW3 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zisolW3
pzisolW3 <- 2*pnorm(-abs(zisolW3))
pzisolW3

# for path 1
# Calculate membership in path  OP*SM*CM 
path1isolW3 <- compute("OP*SM*CM ", data=tff)

obspath1isolW3 <- as.numeric(path1isolW3 <= tff$W)
zpath1isolW3 <- round(((((sum(obspath1isolW3 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zpath1isolW3
pzpath1isolW3 <- 2*pnorm(-abs(zpath1isolW3))
pzpath1isolW3

#for path 2
# Calculate membership in path  TP*SM*CM 
path2isolW3 <- compute("TP*SM*CM ", data=tff)

obspath2isolW3 <- as.numeric(path2isolW3 <= tff$W)
zpath2isolW3 <- round(((((sum(obspath2isolW3 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zpath2isolW3
pzpath2isolW3 <- 2*pnorm(-abs(zpath2isolW3))
pzpath2isolW3

# for path 3
# Calculate membership in path  SP*tp*op*CM
path3isolW3 <- compute("SP*tp*op*CM ", data=tff)

obspath3isolW3 <- as.numeric(path3isolW3 <= tff$W)
zpath3isolW3 <- round(((((sum(obspath3isolW3 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zpath3isolW3
pzpath3isolW3 <- 2*pnorm(-abs(zpath3isolW3))
pzpath3isolW3


####W4 frequency threshold 11, raw consistency threshold 0.929####


ttW4 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.929, n.cut=11, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW4


colnames(tff)

#parsimonious solution
psW4 <- eqmcc(ttW4, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW4

#exclude contradictory assumptions & rows: sp*tp*op (contradict statement of necessity)

ttW4new <- ttW4
nec
rows <- findRows("sp*tp*op", ttW4new, remainders = FALSE)
rows

ttW4new$tt[as.character(rows), "OUT"] <- 0
ttW4new


#conservative solution
csW4 <- eqmcc(ttW4new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW4

#enhanced parsimonious solution (all simplifying assumptions)

epsW4 <- eqmcc(ttW4new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW4


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term TP*CM + (SP*CM)
epsolW4 <- compute("TP*CM + (SP*CM)", data=tff)

obsepsolW4 <- as.numeric(epsolW4 <= tff$W)
zepsolW4 <- round(((((sum(obsepsolW4 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zepsolW4
pzepsolW4 <- 2*pnorm(-abs(zepsolW4))
pzepsolW4



#intermediate solution with directional expectations SM -> W, CM -> W
isW4 <- eqmcc(ttW4new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW4

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
epsW4

# calculate membership in paths of M2: TP*CM + (SP*CM)

path1W4 <- compute("TP*CM ", data=tff)
path2W4 <- compute("SP*CM ", data=tff)


A1path1W4 <- round(sum(as.numeric((path1W4 <= tff$W) & (path1W4 >= nw)))/(1087/100), digits = 1)

A2path1W4 <- round(sum(as.numeric((path1W4 >= tff$W)&(path1W4 >= nw)))/(1087/100), digits = 1)

A3path1W4 <- round(sum(as.numeric((path1W4 >= tff$W)&(path1W4 <= nw)))/(1087/100), digits=1)

A4path1W4 <- round(sum(as.numeric((path1W4 <= tff$W)&(path1W4 <= nw)))/(1087/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W4, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W4)
text(x=0.8, y=0.5, labels = A2path1W4)
text(x=0.5, y=0.2, labels = A3path1W4)
text(x=0.2, y=0.5, labels = A4path1W4)



A1path2W4 <- round(sum(as.numeric((path2W4 <= tff$W) & (path2W4 >= nw)))/(1087/100), digits = 1)

A2path2W4 <- round(sum(as.numeric((path2W4 >= tff$W)&(path2W4 >= nw)))/(1087/100), digits = 1)

A3path2W4 <- round(sum(as.numeric((path2W4 >= tff$W)&(path2W4 <= nw)))/(1087/100), digits=1)

A4path2W4 <- round(sum(as.numeric((path2W4 <= tff$W)&(path2W4 <= nw)))/(1087/100), digits = 1)

plot(tff$W, path2W4, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab='Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W4)
text(x=0.8, y=0.5, labels = A2path2W4)
text(x=0.5, y=0.2, labels = A3path2W4)
text(x=0.2, y=0.5, labels = A4path2W4)

A1epsolW4 <- round(sum(as.numeric((epsolW4 <= tff$W) & (epsolW4 >= nw)))/(1087/100), digits = 1)

A2epsolW4 <- round(sum(as.numeric((epsolW4 >= tff$W)&(epsolW4 >= nw)))/(1087/100), digits = 1)

A3epsolW4 <- round(sum(as.numeric((epsolW4 >= tff$W)&(epsolW4 <= nw)))/(1087/100), digits=1)

A4epsolW4 <- round(sum(as.numeric((epsolW4 <= tff$W)&(epsolW4 <= nw)))/(1087/100), digits = 1)

isW4
plot(tff$W, epsolW4, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW4)
text(x=0.8, y=0.5, labels = A2epsolW4)
text(x=0.5, y=0.2, labels = A3epsolW4)
text(x=0.2, y=0.5, labels = A4epsolW4)


###### W5: frequency threshold 11, raw consistency threshold 0.926######

ttW5 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.926, n.cut=11, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW5


colnames(tff)

#parsimonious solution
psW5 <- eqmcc(ttW5, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW5

#exclude contradictory assumptions & rows: sp*tp*op (contradict statement of necessity)

ttW5new <- ttW5
nec
rows <- findRows("sp*tp*op", ttW5new, remainders = FALSE)
rows

ttW5new$tt[as.character(rows), "OUT"] <- 0
ttW5new


#conservative solution
csW5 <- eqmcc(ttW5new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW5

#enhanced parsimonious solution (all simplifying assumptions)

epsW5 <- eqmcc(ttW5new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW5

#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term SP*CM + sp*TP*SM
epsolW5 <- compute("SP*CM + sp*TP*SM", data=tff)

obsepsolW5 <- as.numeric(epsolW5 <= tff$W)
zepsolW5 <- round(((((sum(obsepsolW5 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zepsolW5
pzepsolW5 <- 2*pnorm(-abs(zepsolW5))
pzepsolW5


#intermediate solution with directional expectations SM -> W, CM -> W
isW5 <- eqmcc(ttW5new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW5



#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: SP*CM + sp*TP*SM

path1W5 <- compute("SP*CM ", data=tff)
path2W5 <- compute("sp*TP*SM ", data=tff)

A1path1W5 <- round(sum(as.numeric((path1W5 <= tff$W) & (path1W5 >= nw)))/(1087/100), digits = 1)

A2path1W5 <- round(sum(as.numeric((path1W5 >= tff$W)&(path1W5 >= nw)))/(1087/100), digits = 1)

A3path1W5 <- round(sum(as.numeric((path1W5 >= tff$W)&(path1W5 <= nw)))/(1087/100), digits=1)

A4path1W5 <- round(sum(as.numeric((path1W5 <= tff$W)&(path1W5 <= nw)))/(1087/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W5)
text(x=0.8, y=0.5, labels = A2path1W5)
text(x=0.5, y=0.2, labels = A3path1W5)
text(x=0.2, y=0.5, labels = A4path1W5)



A1path2W5 <- round(sum(as.numeric((path2W5 <= tff$W) & (path2W5 >= nw)))/(1087/100), digits = 1)

A2path2W5 <- round(sum(as.numeric((path2W5 >= tff$W)&(path2W5 >= nw)))/(1087/100), digits = 1)

A3path2W5 <- round(sum(as.numeric((path2W5 >= tff$W)&(path2W5 <= nw)))/(1087/100), digits=1)

A4path2W5 <- round(sum(as.numeric((path2W5 <= tff$W)&(path2W5 <= nw)))/(1087/100), digits = 1)

plot(tff$W, path2W5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab='Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W5)
text(x=0.8, y=0.5, labels = A2path2W5)
text(x=0.5, y=0.2, labels = A3path2W5)
text(x=0.2, y=0.5, labels = A4path2W5)

A1epsolW5 <- round(sum(as.numeric((epsolW5 <= tff$W) & (epsolW5 >= nw)))/(1087/100), digits = 1)

A2epsolW5 <- round(sum(as.numeric((epsolW5 >= tff$W)&(epsolW5 >= nw)))/(1087/100), digits = 1)

A3epsolW5 <- round(sum(as.numeric((epsolW5 >= tff$W)&(epsolW5 <= nw)))/(1087/100), digits=1)

A4epsolW5 <- round(sum(as.numeric((epsolW5 <= tff$W)&(epsolW5 <= nw)))/(1087/100), digits = 1)

isW5
plot(tff$W, epsolW5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW5)
text(x=0.8, y=0.5, labels = A2epsolW5)
text(x=0.5, y=0.2, labels = A3epsolW5)
text(x=0.2, y=0.5, labels = A4epsolW5)



####W6 frequency threshold 11, raw consistency threshold 0.923####
ttW6 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.923, n.cut=11, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW6


#parsimonious solution
psW6 <- eqmcc(ttW6, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW6

#exclude contradictory assumptions & rows: sp*tp*op  (contradict statement of necessity)

ttW6new <- ttW6
nec
rows <- findRows("sp*tp*op ", ttW6new, remainders = FALSE)
rows

ttW6new$tt[as.character(rows), "OUT"] <- 0
ttW6new


#conservative solution
csW6 <- eqmcc(ttW6new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW6

#enhanced parsimonious solution (all simplifying assumptions)

epsW6 <- eqmcc(ttW6new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW6


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term SP*CM + sp*TP*SM + TP*OP*SM + SP*tp*op*SM
epsolW6 <- compute("SP*CM + sp*TP*SM + TP*OP*SM + SP*tp*op*SM", data=tff)

obsepsolW6 <- as.numeric(epsolW6 <= tff$W)
zepsolW6 <- round(((((sum(obsepsolW6 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zepsolW6
pzepsolW6 <- 2*pnorm(-abs(zepsolW6))
pzepsolW6




#intermediate solution with directional expectations SM -> W, CM -> W
isW6 <- eqmcc(ttW6new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW6

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: SP*CM + sp*TP*SM + TP*OP*SM + SP*tp*op*SM

path1W6 <- compute("SP*CM ", data=tff)
path2W6 <- compute("sp*TP*SM ", data=tff)
path3W6 <- compute("TP*OP*SM ", data=tff)
path4W6 <- compute("SP*tp*op*SM ", data=tff)

A1path1W6 <- round(sum(as.numeric((path1W6 <= tff$W) & (path1W6 >= nw)))/(1087/100), digits = 1)

A2path1W6 <- round(sum(as.numeric((path1W6 >= tff$W)&(path1W6 >= nw)))/(1087/100), digits = 1)

A3path1W6 <- round(sum(as.numeric((path1W6 >= tff$W)&(path1W6 <= nw)))/(1087/100), digits=1)

A4path1W6 <- round(sum(as.numeric((path1W6 <= tff$W)&(path1W6 <= nw)))/(1087/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W6)
text(x=0.8, y=0.5, labels = A2path1W6)
text(x=0.5, y=0.2, labels = A3path1W6)
text(x=0.2, y=0.5, labels = A4path1W6)



A1path2W6 <- round(sum(as.numeric((path2W6 <= tff$W) & (path2W6 >= nw)))/(1087/100), digits = 1)

A2path2W6 <- round(sum(as.numeric((path2W6 >= tff$W)&(path2W6 >= nw)))/(1087/100), digits = 1)

A3path2W6 <- round(sum(as.numeric((path2W6 >= tff$W)&(path2W6 <= nw)))/(1087/100), digits=1)

A4path2W6 <- round(sum(as.numeric((path2W6 <= tff$W)&(path2W6 <= nw)))/(1087/100), digits = 1)

plot(tff$W, path2W6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab='Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W6)
text(x=0.8, y=0.5, labels = A2path2W6)
text(x=0.5, y=0.2, labels = A3path2W6)
text(x=0.2, y=0.5, labels = A4path2W6)

A1path3W6 <- round(sum(as.numeric((path3W6 <= tff$W) & (path3W6 >= nw)))/(1087/100), digits = 1)

A2path3W6 <- round(sum(as.numeric((path3W6 >= tff$W)&(path3W6 >= nw)))/(1087/100), digits = 1)

A3path3W6 <- round(sum(as.numeric((path3W6 >= tff$W)&(path3W6 <= nw)))/(1087/100), digits=1)

A4path3W6 <- round(sum(as.numeric((path3W6 <= tff$W)&(path3W6 <= nw)))/(1087/100), digits = 1)


plot(tff$W, path3W6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 3',
     xlab= 'Path 3',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3W6)
text(x=0.8, y=0.5, labels = A2path3W6)
text(x=0.5, y=0.2, labels = A3path3W6)
text(x=0.2, y=0.5, labels = A4path3W6)


A1path4W6 <- round(sum(as.numeric((path4W6 <= tff$W) & (path4W6 >= nw)))/(1087/100), digits = 1)

A2path4W6 <- round(sum(as.numeric((path4W6 >= tff$W)&(path4W6 >= nw)))/(1087/100), digits = 1)

A3path4W6 <- round(sum(as.numeric((path4W6 >= tff$W)&(path4W6 <= nw)))/(1087/100), digits=1)

A4path4W6 <- round(sum(as.numeric((path4W6 <= tff$W)&(path4W6 <= nw)))/(1087/100), digits = 1)

plot(tff$W, path4W6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 4',
     xlab= 'Path 4',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path4W6)
text(x=0.8, y=0.5, labels = A2path4W6)
text(x=0.5, y=0.2, labels = A3path4W6)
text(x=0.2, y=0.5, labels = A4path4W6)



A1epsolW6 <- round(sum(as.numeric((epsolW6 <= tff$W) & (epsolW6 >= nw)))/(1087/100), digits = 1)

A2epsolW6 <- round(sum(as.numeric((epsolW6 >= tff$W)&(epsolW6 >= nw)))/(1087/100), digits = 1)

A3epsolW6 <- round(sum(as.numeric((epsolW6 >= tff$W)&(epsolW6 <= nw)))/(1087/100), digits=1)

A4epsolW6 <- round(sum(as.numeric((epsolW6 <= tff$W)&(epsolW6 <= nw)))/(1087/100), digits = 1)

isW6
plot(tff$W, epsolW6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW6)
text(x=0.8, y=0.5, labels = A2epsolW6)
text(x=0.5, y=0.2, labels = A3epsolW6)
text(x=0.2, y=0.5, labels = A4epsolW6)


#####W7 frequency threshold 54, raw consistency threshold 0.891####

ttW7 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.891, n.cut=54, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW7


#parsimonious solution
psW7 <- eqmcc(ttW7, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW7

#exclude contradictory assumptions & rows: sp*tp*op  (contradict statement of necessity)

ttW7new <- ttW7
nec
rows <- findRows("sp*tp*op", ttW7new, remainders = FALSE)
rows

ttW7new$tt[as.character(rows), "OUT"] <- 0
ttW7new


#conservative solution
csW7 <- eqmcc(ttW7new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW7

#enhanced parsimonious solution (all simplifying assumptions)

epsW7 <- eqmcc(ttW7new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW7


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term OP*CM
epsolW7 <- compute("OP*CM", data=tff)

obsepsolW7 <- as.numeric(epsolW7 <= tff$W)
zepsolW7 <- round(((((sum(obsepsolW7 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zepsolW7
pzepsolW7 <- 2*pnorm(-abs(zepsolW7))
pzepsolW7

#intermediate solution with directional expectations SM -> W, CM -> W
isW7 <- eqmcc(ttW7new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW7

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: SP*TP*SM

epsW7

path1W7 <- compute("OP*CM", data=tff)

A1path1W7 <- round(sum(as.numeric((path1W7 <= tff$W) & (path1W7 >= nw)))/(1087/100), digits = 1)

A2path1W7 <- round(sum(as.numeric((path1W7 >= tff$W)&(path1W7 >= nw)))/(1087/100), digits = 1)

A3path1W7 <- round(sum(as.numeric((path1W7 >= tff$W)&(path1W7 <= nw)))/(1087/100), digits=1)

A4path1W7 <- round(sum(as.numeric((path1W7 <= tff$W)&(path1W7 <= nw)))/(1087/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W7, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W7)
text(x=0.8, y=0.5, labels = A2path1W7)
text(x=0.5, y=0.2, labels = A3path1W7)
text(x=0.2, y=0.5, labels = A4path1W7)


A1epsolW7 <- round(sum(as.numeric((epsolW7 <= tff$W) & (epsolW7 >= nw)))/(1087/100), digits = 1)

A2epsolW7 <- round(sum(as.numeric((epsolW7 >= tff$W)&(epsolW7 >= nw)))/(1087/100), digits = 1)

A3epsolW7 <- round(sum(as.numeric((epsolW7 >= tff$W)&(epsolW7 <= nw)))/(1087/100), digits=1)

A4epsolW7 <- round(sum(as.numeric((epsolW7 <= tff$W)&(epsolW7 <= nw)))/(1087/100), digits = 1)


plot(tff$W, epsolW7, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW7)
text(x=0.8, y=0.5, labels = A2epsolW7)
text(x=0.5, y=0.2, labels = A3epsolW7)
text(x=0.2, y=0.5, labels = A4epsolW7)


#####W8 frequency threshold 54, raw consistency threshold 0.884####

ttW8 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.884, n.cut=54, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW8


colnames(tff)

#parsimonious solution
psW8 <- eqmcc(ttW8, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW8

#exclude contradictory assumptions & rows: sp*tp*op (contradict statement of necessity)

ttW8new <- ttW8
nec
rows <- findRows("sp*tp*op", ttW8new, remainders = FALSE)
rows

ttW8new$tt[as.character(rows), "OUT"] <- 0
ttW8new


#conservative solution
csW8 <- eqmcc(ttW8new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW8

#enhanced parsimonious solution (all simplifying assumptions)

epsW8 <- eqmcc(ttW8new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW8


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term SP*TP*op + (OP*CM)
epsolW8 <- compute("SP*TP*op + (OP*CM)", data=tff)

obsepsolW8 <- as.numeric(epsolW8 <= tff$W)
zepsolW8 <- round(((((sum(obsepsolW8 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zepsolW8
pzepsolW8 <- 2*pnorm(-abs(zepsolW8))
pzepsolW8

#intermediate solution with directional expectations SM -> W, CM -> W
isW8 <- eqmcc(ttW8new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW8

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: SP*TP*op + (OP*CM)

epsW8

path1W8 <- compute("SP*TP*op ", data=tff)
path2W8 <- compute("OP*CM ", data=tff)


A1path1W8 <- round(sum(as.numeric((path1W8 <= tff$W) & (path1W8 >= nw)))/(1087/100), digits = 1)

A2path1W8 <- round(sum(as.numeric((path1W8 >= tff$W)&(path1W8 >= nw)))/(1087/100), digits = 1)

A3path1W8 <- round(sum(as.numeric((path1W8 >= tff$W)&(path1W8 <= nw)))/(1087/100), digits=1)

A4path1W8 <- round(sum(as.numeric((path1W8 <= tff$W)&(path1W8 <= nw)))/(1087/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W8, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W8)
text(x=0.8, y=0.5, labels = A2path1W8)
text(x=0.5, y=0.2, labels = A3path1W8)
text(x=0.2, y=0.5, labels = A4path1W8)

A1path2W8 <- round(sum(as.numeric((path2W8 <= tff$W) & (path2W8 >= nw)))/(1087/100), digits = 1)

A2path2W8 <- round(sum(as.numeric((path2W8 >= tff$W)&(path2W8 >= nw)))/(1087/100), digits = 1)

A3path2W8 <- round(sum(as.numeric((path2W8 >= tff$W)&(path2W8 <= nw)))/(1087/100), digits=1)

A4path2W8 <- round(sum(as.numeric((path2W8 <= tff$W)&(path2W8 <= nw)))/(1087/100), digits = 1)


plot(tff$W, path2W8, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W8)
text(x=0.8, y=0.5, labels = A2path2W8)
text(x=0.5, y=0.2, labels = A3path2W8)
text(x=0.2, y=0.5, labels = A4path2W8)


A1epsolW8 <- round(sum(as.numeric((epsolW8 <= tff$W) & (epsolW8 >= nw)))/(1087/100), digits = 1)

A2epsolW8 <- round(sum(as.numeric((epsolW8 >= tff$W)&(epsolW8 >= nw)))/(1087/100), digits = 1)

A3epsolW8 <- round(sum(as.numeric((epsolW8 >= tff$W)&(epsolW8 <= nw)))/(1087/100), digits=1)

A4epsolW8 <- round(sum(as.numeric((epsolW8 <= tff$W)&(epsolW8 <= nw)))/(1087/100), digits = 1)

isW8
plot(tff$W, epsolW8, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW8)
text(x=0.8, y=0.5, labels = A2epsolW8)
text(x=0.5, y=0.2, labels = A3epsolW8)
text(x=0.2, y=0.5, labels = A4epsolW8)


#####W9 frequency threshold 54, raw consistency threshold 0.883####

ttW9 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.883, n.cut=54, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW9


#parsimonious solution
psW9 <- eqmcc(ttW9, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW9

#exclude contradictory assumptions & rows: sp*tp*op (contradict statement of necessity)

ttW9new <- ttW9
nec
rows <- findRows("sp*tp*op ", ttW9new, remainders = FALSE)
rows

ttW9new$tt[as.character(rows), "OUT"] <- 0
ttW9new


#conservative solution
csW9 <- eqmcc(ttW9new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW9

#enhanced parsimonious solution (all simplifying assumptions)

epsW9 <- eqmcc(ttW9new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW9

#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term sp*OP + SP*TP*op + (SP*CM)
epsolW9 <- compute("sp*OP + SP*TP*op + (SP*CM)", data=tff)

obsepsolW9 <- as.numeric(epsolW9 <= tff$W)
zepsolW9 <- round(((((sum(obsepsolW9 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zepsolW9
pzepsolW9 <- 2*pnorm(-abs(zepsolW9))
pzepsolW9


#intermediate solution with directional expectations SM -> W, CM -> W
isW9 <- eqmcc(ttW9new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW9

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M3: sp*OP + SP*TP*op + (SP*CM)

path1W9 <- compute("sp*OP ", data=tff)
path2W9 <- compute("SP*TP*op ", data=tff)
path3W9 <- compute("SP*CM ", data=tff)

A1path1W9 <- round(sum(as.numeric((path1W9 <= tff$W) & (path1W9 >= nw)))/(1087/100), digits = 1)

A2path1W9 <- round(sum(as.numeric((path1W9 >= tff$W)&(path1W9 >= nw)))/(1087/100), digits = 1)

A3path1W9 <- round(sum(as.numeric((path1W9 >= tff$W)&(path1W9 <= nw)))/(1087/100), digits=1)

A4path1W9 <- round(sum(as.numeric((path1W9 <= tff$W)&(path1W9 <= nw)))/(1087/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W9, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W9)
text(x=0.8, y=0.5, labels = A2path1W9)
text(x=0.5, y=0.2, labels = A3path1W9)
text(x=0.2, y=0.5, labels = A4path1W9)


A1path2W9 <- round(sum(as.numeric((path2W9 <= tff$W) & (path2W9 >= nw)))/(1087/100), digits = 1)

A2path2W9 <- round(sum(as.numeric((path2W9 >= tff$W)&(path2W9 >= nw)))/(1087/100), digits = 1)

A3path2W9 <- round(sum(as.numeric((path2W9 >= tff$W)&(path2W9 <= nw)))/(1087/100), digits=1)

A4path2W9 <- round(sum(as.numeric((path2W9 <= tff$W)&(path2W9 <= nw)))/(1087/100), digits = 1)


plot(tff$W, path2W9, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W9)
text(x=0.8, y=0.5, labels = A2path2W9)
text(x=0.5, y=0.2, labels = A3path2W9)
text(x=0.2, y=0.5, labels = A4path2W9)

A1path3W9 <- round(sum(as.numeric((path3W9 <= tff$W) & (path3W9 >= nw)))/(1087/100), digits = 1)

A2path3W9 <- round(sum(as.numeric((path3W9 >= tff$W)&(path3W9 >= nw)))/(1087/100), digits = 1)

A3path3W9 <- round(sum(as.numeric((path3W9 >= tff$W)&(path3W9 <= nw)))/(1087/100), digits=1)

A4path3W9 <- round(sum(as.numeric((path3W9 <= tff$W)&(path3W9 <= nw)))/(1087/100), digits = 1)


plot(tff$W, path3W9, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 3',
     xlab= 'Path 3',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3W9)
text(x=0.8, y=0.5, labels = A2path3W9)
text(x=0.5, y=0.2, labels = A3path3W9)
text(x=0.2, y=0.5, labels = A4path3W9)




A1epsolW9 <- round(sum(as.numeric((epsolW9 <= tff$W) & (epsolW9 >= nw)))/(1087/100), digits = 1)

A2epsolW9 <- round(sum(as.numeric((epsolW9 >= tff$W)&(epsolW9 >= nw)))/(1087/100), digits = 1)

A3epsolW9 <- round(sum(as.numeric((epsolW9 >= tff$W)&(epsolW9 <= nw)))/(1087/100), digits=1)

A4epsolW9 <- round(sum(as.numeric((epsolW9 <= tff$W)&(epsolW9 <= nw)))/(1087/100), digits = 1)

isW9
plot(tff$W, epsolW9, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW9)
text(x=0.8, y=0.5, labels = A2epsolW9)
text(x=0.5, y=0.2, labels = A3epsolW9)
text(x=0.2, y=0.5, labels = A4epsolW9)



#######################ANALYSIS OF SUFFICIENCY FOR w ####################


tff <- read.csv("tff_d2_c1.csv", row.names=1, header = TRUE, sep = ";",  dec = ",")
describe(tff)

nw <- 1-tff$W

# unenable assumptions: 

suf <- "OP*SM*CM + SP*SM*CM + TP*SM*CM"
suf

#####w1  frequency threshold 5, raw consistency threshold 0.882 #####


ttw1 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.882, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw1


#parsimonious solution
psw1 <- eqmcc(ttw1, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw1

#exclude contradictory assumptions & rows: OP*SM*CM + SP*SM*CM + TP*SM*CM (contradict statement of sufficiency)

ttw1new <- ttw1
suf
rows <- findRows("OP*SM*CM + SP*SM*CM + TP*SM*CM", ttw1new, remainders = FALSE)
rows

ttw1new$tt[as.character(rows), "OUT"] <- 0
ttw1new



#conservative solution
csw1 <- eqmcc(ttw1new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw1

#enhanced parsimonious solution (all simplifying assumptions)

epsw1 <- eqmcc(ttw1new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw1


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  op*sm*CM + sp*sm*CM + sp*tp*op*CM + sp*tp*OP*SM*cm
epsolw1 <- compute("op*sm*CM + sp*sm*CM + sp*tp*op*CM + sp*tp*OP*SM*cm ", data=tff)

obsepsolw1 <- as.numeric(epsolw1 <= nw)
zepsolw1 <- round(((((sum(obsepsolw1 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zepsolw1
pzepsolw1 <- 2*pnorm(-abs(zepsolw1))
pzepsolw1



#intermediate solution with Directional expectations: sp->w, tp->w, op->w.

isw1 <- eqmcc(ttw1new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw1
epsw1


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: tp*op*CM + sp*op*sm*CM
path1w1 <- compute("op*sm*CM ", data=tff)
path2w1 <- compute("sp*sm*CM ", data=tff)
path3w1 <- compute("sp*tp*op*CM ", data=tff)
path4w1 <- compute("sp*tp*OP*SM*cm", data=tff)

A1path1w1 <- round(sum(as.numeric((path1w1 <= nw) & (path1w1 >= tff$W)))/(1087/100), digits = 1)

A2path1w1 <- round(sum(as.numeric((path1w1 >= nw)&(path1w1 >= tff$W)))/(1087/100), digits = 1)

A3path1w1 <- round(sum(as.numeric((path1w1 >= nw)&(path1w1 <= tff$W)))/(1087/100), digits=1)

A4path1w1 <- round(sum(as.numeric((path1w1 <= nw)&(path1w1 <= tff$W)))/(1087/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w1)
text(x=0.8, y=0.5, labels = A2path1w1)
text(x=0.5, y=0.2, labels = A3path1w1)
text(x=0.2, y=0.5, labels = A4path1w1)



A1path2w1 <- round(sum(as.numeric((path2w1 <= nw) & (path2w1 >= tff$W)))/(1087/100), digits = 1)

A2path2w1 <- round(sum(as.numeric((path2w1 >= nw)&(path2w1 >= tff$W)))/(1087/100), digits = 1)

A3path2w1 <- round(sum(as.numeric((path2w1 >= nw)&(path2w1 <= tff$W)))/(1087/100), digits=1)

A4path2w1 <- round(sum(as.numeric((path2w1 <= nw)&(path2w1 <= tff$W)))/(1087/100), digits = 1)


plot(nw, path2w1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w1)
text(x=0.8, y=0.5, labels = A2path2w1)
text(x=0.5, y=0.2, labels = A3path2w1)
text(x=0.2, y=0.5, labels = A4path2w1)


A1path3w1 <- round(sum(as.numeric((path3w1 <= nw) & (path3w1 >= tff$W)))/(1087/100), digits = 1)

A2path3w1 <- round(sum(as.numeric((path3w1 >= nw)&(path3w1 >= tff$W)))/(1087/100), digits = 1)

A3path3w1 <- round(sum(as.numeric((path3w1 >= nw)&(path3w1 <= tff$W)))/(1087/100), digits=1)

A4path3w1 <- round(sum(as.numeric((path3w1 <= nw)&(path3w1 <= tff$W)))/(1087/100), digits = 1)


plot(nw, path3w1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 3',
     xlab= 'Path 3',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3w1)
text(x=0.8, y=0.5, labels = A2path3w1)
text(x=0.5, y=0.2, labels = A3path3w1)
text(x=0.2, y=0.5, labels = A4path3w1)


A1path4w1 <- round(sum(as.numeric((path4w1 <= nw) & (path4w1 >= tff$W)))/(1087/100), digits = 1)

A2path4w1 <- round(sum(as.numeric((path4w1 >= nw)&(path4w1 >= tff$W)))/(1087/100), digits = 1)

A3path4w1 <- round(sum(as.numeric((path4w1 >= nw)&(path4w1 <= tff$W)))/(1087/100), digits=1)

A4path4w1 <- round(sum(as.numeric((path4w1 <= nw)&(path4w1 <= tff$W)))/(1087/100), digits = 1)


plot(nw, path4w1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 4',
     xlab= 'Path 4',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path4w1)
text(x=0.8, y=0.5, labels = A2path4w1)
text(x=0.5, y=0.2, labels = A3path4w1)
text(x=0.2, y=0.5, labels = A4path4w1)


A1epsolw1 <- round(sum(as.numeric((epsolw1 <= nw) & (epsolw1 >= tff$W)))/(1087/100), digits = 1)

A2epsolw1 <- round(sum(as.numeric((epsolw1 >= nw)&(epsolw1 >= tff$W)))/(1087/100), digits = 1)

A3epsolw1 <- round(sum(as.numeric((epsolw1 >= nw)&(epsolw1 <= tff$W)))/(1087/100), digits=1)

A4epsolw1 <- round(sum(as.numeric((epsolw1 <= nw)&(epsolw1 <= tff$W)))/(1087/100), digits = 1)


plot(nw, epsolw1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw1)
text(x=0.8, y=0.5, labels = A2epsolw1)
text(x=0.5, y=0.2, labels = A3epsolw1)
text(x=0.2, y=0.5, labels = A4epsolw1)



#####w2  frequency threshold 5, raw consistency threshold 0.880 #####


ttw2 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.880, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw2

#parsimonious solution
psw2 <- eqmcc(ttw2, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw2

#exclude contradictory assumptions & rows: OP*SM*CM + SP*SM*CM + TP*SM*CM (contradict statement of sufficiency)

ttw2new <- ttw2
suf
rows <- findRows("OP*SM*CM + SP*SM*CM + TP*SM*CM", ttw2new, remainders = FALSE)
rows

ttw2new$tt[as.character(rows), "OUT"] <- 0
ttw2new


#conservative solution
csw2 <- eqmcc(ttw2new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw2

#enhanced parsimonious solution (all simplifying assumptions)

epsw2 <- eqmcc(ttw2new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw2

#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  op*sm*CM + sp*sm*CM + sp*tp*op*CM + sp*tp*OP*SM*cm + SP*tp*op*SM*cm
epsolw2 <- compute("op*sm*CM + sp*sm*CM + sp*tp*op*CM + sp*tp*OP*SM*cm + SP*tp*op*SM*cm ", data=tff)

obsepsolw2 <- as.numeric(epsolw2 <= nw)
zepsolw2 <- round(((((sum(obsepsolw2 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zepsolw2
pzepsolw2 <- 2*pnorm(-abs(zepsolw2))
pzepsolw2


#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw2 <- eqmcc(ttw2new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw2
epsw2

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: op*sm*CM + sp*sm*CM + sp*tp*op*CM + sp*tp*OP*SM*cm + SP*tp*op*SM*cm



path1w2 <- compute("op*sm*CM ", data=tff)
path2w2 <- compute("sp*sm*CM ", data=tff)
path3w2 <- compute("sp*tp*op*CM ", data=tff)
path4w2 <- compute("sp*tp*OP*SM*cm", data=tff)
path5w2 <- compute("SP*tp*op*SM*cm", data=tff)

A1path1w2 <- round(sum(as.numeric((path1w2 <= nw) & (path1w2 >= tff$W)))/(1087/100), digits = 1)

A2path1w2 <- round(sum(as.numeric((path1w2 >= nw)&(path1w2 >= tff$W)))/(1087/100), digits = 1)

A3path1w2 <- round(sum(as.numeric((path1w2 >= nw)&(path1w2 <= tff$W)))/(1087/100), digits=1)

A4path1w2 <- round(sum(as.numeric((path1w2 <= nw)&(path1w2 <= tff$W)))/(1087/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w2)
text(x=0.8, y=0.5, labels = A2path1w2)
text(x=0.5, y=0.2, labels = A3path1w2)
text(x=0.2, y=0.5, labels = A4path1w2)



A1path2w2 <- round(sum(as.numeric((path2w2 <= nw) & (path2w2 >= tff$W)))/(1087/100), digits = 1)

A2path2w2 <- round(sum(as.numeric((path2w2 >= nw)&(path2w2 >= tff$W)))/(1087/100), digits = 1)

A3path2w2 <- round(sum(as.numeric((path2w2 >= nw)&(path2w2 <= tff$W)))/(1087/100), digits=1)

A4path2w2 <- round(sum(as.numeric((path2w2 <= nw)&(path2w2 <= tff$W)))/(1087/100), digits = 1)


plot(nw, path2w2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w2)
text(x=0.8, y=0.5, labels = A2path2w2)
text(x=0.5, y=0.2, labels = A3path2w2)
text(x=0.2, y=0.5, labels = A4path2w2)


A1path3w2 <- round(sum(as.numeric((path3w2 <= nw) & (path3w2 >= tff$W)))/(1087/100), digits = 1)

A2path3w2 <- round(sum(as.numeric((path3w2 >= nw)&(path3w2 >= tff$W)))/(1087/100), digits = 1)

A3path3w2 <- round(sum(as.numeric((path3w2 >= nw)&(path3w2 <= tff$W)))/(1087/100), digits=1)

A4path3w2 <- round(sum(as.numeric((path3w2 <= nw)&(path3w2 <= tff$W)))/(1087/100), digits = 1)


plot(nw, path3w2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 3',
     xlab= 'Path 3',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3w2)
text(x=0.8, y=0.5, labels = A2path3w2)
text(x=0.5, y=0.2, labels = A3path3w2)
text(x=0.2, y=0.5, labels = A4path3w2)


A1path4w2 <- round(sum(as.numeric((path4w2 <= nw) & (path4w2 >= tff$W)))/(1087/100), digits = 1)

A2path4w2 <- round(sum(as.numeric((path4w2 >= nw)&(path4w2 >= tff$W)))/(1087/100), digits = 1)

A3path4w2 <- round(sum(as.numeric((path4w2 >= nw)&(path4w2 <= tff$W)))/(1087/100), digits=1)

A4path4w2 <- round(sum(as.numeric((path4w2 <= nw)&(path4w2 <= tff$W)))/(1087/100), digits = 1)


plot(nw, path4w2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 4',
     xlab= 'Path 4',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path4w2)
text(x=0.8, y=0.5, labels = A2path4w2)
text(x=0.5, y=0.2, labels = A3path4w2)
text(x=0.2, y=0.5, labels = A4path4w2)

A1path5w2 <- round(sum(as.numeric((path5w2 <= nw) & (path5w2 >= tff$W)))/(1087/100), digits = 1)

A2path5w2 <- round(sum(as.numeric((path5w2 >= nw)&(path5w2 >= tff$W)))/(1087/100), digits = 1)

A3path5w2 <- round(sum(as.numeric((path5w2 >= nw)&(path5w2 <= tff$W)))/(1087/100), digits=1)

A4path5w2 <- round(sum(as.numeric((path5w2 <= nw)&(path5w2 <= tff$W)))/(1087/100), digits = 1)


plot(nw, path5w2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 5',
     xlab= 'Path 5',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path5w2)
text(x=0.8, y=0.5, labels = A2path5w2)
text(x=0.5, y=0.2, labels = A3path5w2)
text(x=0.2, y=0.5, labels = A4path5w2)


A1epsolw2 <- round(sum(as.numeric((epsolw2 <= nw) & (epsolw2 >= tff$W)))/(1087/100), digits = 1)

A2epsolw2 <- round(sum(as.numeric((epsolw2 >= nw)&(epsolw2 >= tff$W)))/(1087/100), digits = 1)

A3epsolw2 <- round(sum(as.numeric((epsolw2 >= nw)&(epsolw2 <= tff$W)))/(1087/100), digits=1)

A4epsolw2 <- round(sum(as.numeric((epsolw2 <= nw)&(epsolw2 <= tff$W)))/(1087/100), digits = 1)


plot(nw, epsolw2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw2)
text(x=0.8, y=0.5, labels = A2epsolw2)
text(x=0.5, y=0.2, labels = A3epsolw2)
text(x=0.2, y=0.5, labels = A4epsolw2)

#Z-Test for intermediate solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  op*sm*CM + sp*sm*CM + sp*tp*op*CM + sp*tp*OP*SM*cm + SP*tp*op*SM*cm

isolw2 <- compute("op*sm*CM + sp*sm*CM + sp*tp*op*CM + sp*tp*OP*SM*cm + SP*tp*op*SM*cm ", data=tff)

obsisolw2 <- as.numeric(isolw2 <= nw)
zisolw2 <- round(((((sum(obsisolw2 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zisolw2
pzisolw2 <- 2*pnorm(-abs(zisolw2))
pzisolw2

# for path 1
# Calculate membership in path  op*sm*CM 
path1isolw2 <- compute("op*sm*CM ", data=tff)

obspath1isolw2 <- as.numeric(path1isolw2 <= nw)
zpath1isolw2 <- round(((((sum(obspath1isolw2 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zpath1isolw2
pzpath1isolw2 <- 2*pnorm(-abs(zpath1isolw2))
pzpath1isolw2

#for path 2
# Calculate membership in path  sp*sm*CM 
path2isolw2 <- compute("sp*sm*CM ", data=tff)

obspath2isolw2 <- as.numeric(path2isolw2 <= nw)
zpath2isolw2 <- round(((((sum(obspath2isolw2 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zpath2isolw2
pzpath2isolw2 <- 2*pnorm(-abs(zpath2isolw2))
pzpath2isolw2

# for path 3
# Calculate membership in path  sp*tp*op*CM 

path3isolw2 <- compute("sp*tp*op*CM ", data=tff)

obspath3isolw2 <- as.numeric(path3isolw2 <= nw)
zpath3isolw2 <- round(((((sum(obspath3isolw2 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zpath3isolw2
pzpath3isolw2 <- 2*pnorm(-abs(zpath3isolw2))
pzpath3isolw2


# for path 4
# Calculate membership in path  sp*tp*OP*SM*cm 


path4isolw2 <- compute("sp*tp*OP*SM*cm ", data=tff)

obspath4isolw2 <- as.numeric(path4isolw2 <= nw)
zpath4isolw2 <- round(((((sum(obspath4isolw2 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zpath4isolw2
pzpath4isolw2 <- 2*pnorm(-abs(zpath4isolw2))
pzpath4isolw2

# for path 5 SP*tp*op*SM*cm
# Calculate membership in path  SP*tp*op*SM*cm


path5isolw2 <- compute("SP*tp*op*SM*cm ", data=tff)

obspath5isolw2 <- as.numeric(path5isolw2 <= nw)
zpath5isolw2 <- round(((((sum(obspath5isolw2 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zpath5isolw2
pzpath5isolw2 <- 2*pnorm(-abs(zpath5isolw2))
pzpath5isolw2



# identify simplifying assumptions and easy counterfactuals

SAw2 <- psw2$SA
SAw2

ECw2 <- isw2$i.sol$C1P1$EC
ECw2

#####w3  frequency threshold 5, raw consistency threshold 0.8785 #####


ttw3 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.8785, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw3

#parsimonious solution
psw3 <- eqmcc(ttw3, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw3

#exclude contradictory assumptions & rows: OP*SM*CM + SP*SM*CM + TP*SM*CM (contradict statement of sufficiency)

ttw3new <- ttw3
suf
rows <- findRows("OP*SM*CM + SP*SM*CM + TP*SM*CM", ttw3new, remainders = FALSE)
rows

ttw3new$tt[as.character(rows), "OUT"] <- 0
ttw3new


#conservative solution
csw3 <- eqmcc(ttw3new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw3

#enhanced parsimonious solution (all simplifying assumptions)

epsw3 <- eqmcc(ttw3new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw3

#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  op*sm*CM + sp*sm*CM + sp*tp*op*CM + sp*tp*OP*SM*cm + SP*tp*op*SM*cm
epsolw3 <- compute("op*sm*CM + sp*sm*CM + sp*tp*op*CM + sp*tp*OP*SM*cm + SP*tp*op*SM*cm ", data=tff)

obsepsolw3 <- as.numeric(epsolw3 <= nw)
zepsolw3 <- round(((((sum(obsepsolw3 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zepsolw3
pzepsolw3 <- 2*pnorm(-abs(zepsolw3))
pzepsolw3


#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw3 <- eqmcc(ttw3new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw3
epsw3

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: op*sm*CM + sp*sm*CM + sp*tp*op*CM + sp*tp*OP*SM*cm + SP*tp*op*SM*cm



path1w3 <- compute("op*sm*CM ", data=tff)
path2w3 <- compute("sp*sm*CM ", data=tff)
path3w3 <- compute("sp*tp*op*CM ", data=tff)
path4w3 <- compute("sp*tp*OP*SM*cm", data=tff)
path5w3 <- compute("SP*tp*op*SM*cm", data=tff)

A1path1w3 <- round(sum(as.numeric((path1w3 <= nw) & (path1w3 >= tff$W)))/(1087/100), digits = 1)

A2path1w3 <- round(sum(as.numeric((path1w3 >= nw)&(path1w3 >= tff$W)))/(1087/100), digits = 1)

A3path1w3 <- round(sum(as.numeric((path1w3 >= nw)&(path1w3 <= tff$W)))/(1087/100), digits=1)

A4path1w3 <- round(sum(as.numeric((path1w3 <= nw)&(path1w3 <= tff$W)))/(1087/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w3)
text(x=0.8, y=0.5, labels = A2path1w3)
text(x=0.5, y=0.2, labels = A3path1w3)
text(x=0.2, y=0.5, labels = A4path1w3)



A1path2w3 <- round(sum(as.numeric((path2w3 <= nw) & (path2w3 >= tff$W)))/(1087/100), digits = 1)

A2path2w3 <- round(sum(as.numeric((path2w3 >= nw)&(path2w3 >= tff$W)))/(1087/100), digits = 1)

A3path2w3 <- round(sum(as.numeric((path2w3 >= nw)&(path2w3 <= tff$W)))/(1087/100), digits=1)

A4path2w3 <- round(sum(as.numeric((path2w3 <= nw)&(path2w3 <= tff$W)))/(1087/100), digits = 1)


plot(nw, path2w3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w3)
text(x=0.8, y=0.5, labels = A2path2w3)
text(x=0.5, y=0.2, labels = A3path2w3)
text(x=0.2, y=0.5, labels = A4path2w3)


A1path3w3 <- round(sum(as.numeric((path3w3 <= nw) & (path3w3 >= tff$W)))/(1087/100), digits = 1)

A2path3w3 <- round(sum(as.numeric((path3w3 >= nw)&(path3w3 >= tff$W)))/(1087/100), digits = 1)

A3path3w3 <- round(sum(as.numeric((path3w3 >= nw)&(path3w3 <= tff$W)))/(1087/100), digits=1)

A4path3w3 <- round(sum(as.numeric((path3w3 <= nw)&(path3w3 <= tff$W)))/(1087/100), digits = 1)


plot(nw, path3w3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 3',
     xlab= 'Path 3',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3w3)
text(x=0.8, y=0.5, labels = A2path3w3)
text(x=0.5, y=0.2, labels = A3path3w3)
text(x=0.2, y=0.5, labels = A4path3w3)


A1path4w3 <- round(sum(as.numeric((path4w3 <= nw) & (path4w3 >= tff$W)))/(1087/100), digits = 1)

A2path4w3 <- round(sum(as.numeric((path4w3 >= nw)&(path4w3 >= tff$W)))/(1087/100), digits = 1)

A3path4w3 <- round(sum(as.numeric((path4w3 >= nw)&(path4w3 <= tff$W)))/(1087/100), digits=1)

A4path4w3 <- round(sum(as.numeric((path4w3 <= nw)&(path4w3 <= tff$W)))/(1087/100), digits = 1)


plot(nw, path4w3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 4',
     xlab= 'Path 4',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path4w3)
text(x=0.8, y=0.5, labels = A2path4w3)
text(x=0.5, y=0.2, labels = A3path4w3)
text(x=0.2, y=0.5, labels = A4path4w3)

A1path5w3 <- round(sum(as.numeric((path5w3 <= nw) & (path5w3 >= tff$W)))/(1087/100), digits = 1)

A2path5w3 <- round(sum(as.numeric((path5w3 >= nw)&(path5w3 >= tff$W)))/(1087/100), digits = 1)

A3path5w3 <- round(sum(as.numeric((path5w3 >= nw)&(path5w3 <= tff$W)))/(1087/100), digits=1)

A4path5w3 <- round(sum(as.numeric((path5w3 <= nw)&(path5w3 <= tff$W)))/(1087/100), digits = 1)


plot(nw, path5w3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 5',
     xlab= 'Path 5',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path5w3)
text(x=0.8, y=0.5, labels = A2path5w3)
text(x=0.5, y=0.2, labels = A3path5w3)
text(x=0.2, y=0.5, labels = A4path5w3)


A1epsolw3 <- round(sum(as.numeric((epsolw3 <= nw) & (epsolw3 >= tff$W)))/(1087/100), digits = 1)

A2epsolw3 <- round(sum(as.numeric((epsolw3 >= nw)&(epsolw3 >= tff$W)))/(1087/100), digits = 1)

A3epsolw3 <- round(sum(as.numeric((epsolw3 >= nw)&(epsolw3 <= tff$W)))/(1087/100), digits=1)

A4epsolw3 <- round(sum(as.numeric((epsolw3 <= nw)&(epsolw3 <= tff$W)))/(1087/100), digits = 1)


plot(nw, epsolw3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw3)
text(x=0.8, y=0.5, labels = A2epsolw3)
text(x=0.5, y=0.2, labels = A3epsolw3)
text(x=0.2, y=0.5, labels = A4epsolw3)



#####w4  frequency threshold 11, raw consistency threshold 0.880 #####


ttw4 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.880, n.cut=11, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw4

#parsimonious solution
psw4 <- eqmcc(ttw4, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw4

#exclude contradictory assumptions & rows: OP*SM*CM + SP*SM*CM + TP*SM*CM (contradict statement of sufficiency)

ttw4new <- ttw4
suf
rows <- findRows("OP*SM*CM + SP*SM*CM + TP*SM*CM", ttw4new, remainders = FALSE)
rows

ttw4new$tt[as.character(rows), "OUT"] <- 0
ttw4new


#conservative solution
csw4 <- eqmcc(ttw4new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw4

#enhanced parsimonious solution (all simplifying assumptions)

epsw4 <- eqmcc(ttw4new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw4


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  SP*tp*op*SM*cm
epsolw4 <- compute("SP*tp*op*SM*cm ", data=tff)

obsepsolw4 <- as.numeric(epsolw4 <= nw)
zepsolw4 <- round(((((sum(obsepsolw4 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zepsolw4
pzepsolw4 <- 2*pnorm(-abs(zepsolw4))
pzepsolw4

#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw4 <- eqmcc(ttw4new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw4

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: SP*tp*op*SM*cm
path1w4 <- compute("SP*tp*op*SM*cm ", data=tff)

A1path1w4 <- round(sum(as.numeric((path1w4 <= nw) & (path1w4 >= tff$W)))/(1087/100), digits = 1)

A2path1w4 <- round(sum(as.numeric((path1w4 >= nw)&(path1w4 >= tff$W)))/(1087/100), digits = 1)

A3path1w4 <- round(sum(as.numeric((path1w4 >= nw)&(path1w4 <= tff$W)))/(1087/100), digits=1)

A4path1w4 <- round(sum(as.numeric((path1w4 <= nw)&(path1w4 <= tff$W)))/(1087/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w4, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w4)
text(x=0.8, y=0.5, labels = A2path1w4)
text(x=0.5, y=0.2, labels = A3path1w4)
text(x=0.2, y=0.5, labels = A4path1w4)

A1epsolw4 <- round(sum(as.numeric((epsolw4 <= nw) & (epsolw4 >= tff$W)))/(1087/100), digits = 1)

A2epsolw4 <- round(sum(as.numeric((epsolw4 >= nw)&(epsolw4 >= tff$W)))/(1087/100), digits = 1)

A3epsolw4 <- round(sum(as.numeric((epsolw4 >= nw)&(epsolw4 <= tff$W)))/(1087/100), digits=1)

A4epsolw4 <- round(sum(as.numeric((epsolw4 <= nw)&(epsolw4 <= tff$W)))/(1087/100), digits = 1)


plot(nw, epsolw4, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw4)
text(x=0.8, y=0.5, labels = A2epsolw4)
text(x=0.5, y=0.2, labels = A3epsolw4)
text(x=0.2, y=0.5, labels = A4epsolw4)


#####w5  frequency threshold 11, raw consistency threshold 0.868 #####


ttw5 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.868, n.cut=11, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw5

#parsimonious solution
psw5 <- eqmcc(ttw5, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw5

#exclude contradictory assumptions & rows: OP*SM*CM + SP*SM*CM + TP*SM*CM (contradict statement of sufficiency)

ttw5new <- ttw5
suf
rows <- findRows("OP*SM*CM + SP*SM*CM + TP*SM*CM", ttw5new, remainders = FALSE)
rows

ttw5new$tt[as.character(rows), "OUT"] <- 0
ttw5new


#conservative solution
csw5 <- eqmcc(ttw5new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw5

#enhanced parsimonious solution (all simplifying assumptions)

epsw5 <- eqmcc(ttw5new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw5



#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  tp*SM*cm + sp*op*SM*cm
epsolw5 <- compute("tp*SM*cm + sp*op*SM*cm ", data=tff)

obsepsolw5 <- as.numeric(epsolw5 <= nw)
zepsolw5 <- round(((((sum(obsepsolw5 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zepsolw5
pzepsolw5 <- 2*pnorm(-abs(zepsolw5))
pzepsolw5

#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw5 <- eqmcc(ttw5new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw5


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: tp*SM*cm + sp*op*SM*cm

path1w5 <- compute("tp*SM*cm ", data=tff)
path2w5 <- compute("sp*op*SM*cm ", data=tff)

A1path1w5 <- round(sum(as.numeric((path1w5 <= nw) & (path1w5 >= tff$W)))/(1087/100), digits = 1)

A2path1w5 <- round(sum(as.numeric((path1w5 >= nw)&(path1w5 >= tff$W)))/(1087/100), digits = 1)

A3path1w5 <- round(sum(as.numeric((path1w5 >= nw)&(path1w5 <= tff$W)))/(1087/100), digits=1)

A4path1w5 <- round(sum(as.numeric((path1w5 <= nw)&(path1w5 <= tff$W)))/(1087/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w5)
text(x=0.8, y=0.5, labels = A2path1w5)
text(x=0.5, y=0.2, labels = A3path1w5)
text(x=0.2, y=0.5, labels = A4path1w5)

A1path2w5 <- round(sum(as.numeric((path2w5 <= nw) & (path2w5 >= tff$W)))/(1087/100), digits = 1)

A2path2w5 <- round(sum(as.numeric((path2w5 >= nw)&(path2w5 >= tff$W)))/(1087/100), digits = 1)

A3path2w5 <- round(sum(as.numeric((path2w5 >= nw)&(path2w5 <= tff$W)))/(1087/100), digits=1)

A4path2w5 <- round(sum(as.numeric((path2w5 <= nw)&(path2w5 <= tff$W)))/(1087/100), digits = 1)

plot(nw, path2w5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w5)
text(x=0.8, y=0.5, labels = A2path2w5)
text(x=0.5, y=0.2, labels = A3path2w5)
text(x=0.2, y=0.5, labels = A4path2w5)

A1epsolw5 <- round(sum(as.numeric((epsolw5 <= nw) & (epsolw5 >= tff$W)))/(1087/100), digits = 1)

A2epsolw5 <- round(sum(as.numeric((epsolw5 >= nw)&(epsolw5 >= tff$W)))/(1087/100), digits = 1)

A3epsolw5 <- round(sum(as.numeric((epsolw5 >= nw)&(epsolw5 <= tff$W)))/(1087/100), digits=1)

A4epsolw5 <- round(sum(as.numeric((epsolw5 <= nw)&(epsolw5 <= tff$W)))/(1087/100), digits = 1)


plot(nw, epsolw5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw5)
text(x=0.8, y=0.5, labels = A2epsolw5)
text(x=0.5, y=0.2, labels = A3epsolw5)
text(x=0.2, y=0.5, labels = A4epsolw5)


#####w6  frequency threshold 11, raw consistency threshold 0.861#####


ttw6 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.861, n.cut=11, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw6

#parsimonious solution
psw6 <- eqmcc(ttw6, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw6

#exclude contradictory assumptions & rows: OP*SM*CM + SP*SM*CM + TP*SM*CM (contradict statement of sufficiency)

ttw6new <- ttw6
suf
rows <- findRows("OP*SM*CM + SP*SM*CM + TP*SM*CM", ttw6new, remainders = FALSE)
rows

ttw6new$tt[as.character(rows), "OUT"] <- 0
ttw6new


#conservative solution
csw6 <- eqmcc(ttw6new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw6

#enhanced parsimonious solution (all simplifying assumptions)

epsw6 <- eqmcc(ttw6new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw6



#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  tp*SM*cm + sp*op*SM*cm
epsolw6 <- compute("tp*SM*cm + sp*op*SM*cm ", data=tff)

obsepsolw6 <- as.numeric(epsolw6 <= nw)
zepsolw6 <- round(((((sum(obsepsolw6 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zepsolw6
pzepsolw6 <- 2*pnorm(-abs(zepsolw6))
pzepsolw6

#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw6 <- eqmcc(ttw6new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw6


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: tp*SM*cm + sp*op*SM*cm

path1w6 <- compute("tp*SM*cm ", data=tff)
path2w6 <- compute("sp*op*SM*cm ", data=tff)

A1path1w6 <- round(sum(as.numeric((path1w6 <= nw) & (path1w6 >= tff$W)))/(1087/100), digits = 1)

A2path1w6 <- round(sum(as.numeric((path1w6 >= nw)&(path1w6 >= tff$W)))/(1087/100), digits = 1)

A3path1w6 <- round(sum(as.numeric((path1w6 >= nw)&(path1w6 <= tff$W)))/(1087/100), digits=1)

A4path1w6 <- round(sum(as.numeric((path1w6 <= nw)&(path1w6 <= tff$W)))/(1087/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w6)
text(x=0.8, y=0.5, labels = A2path1w6)
text(x=0.5, y=0.2, labels = A3path1w6)
text(x=0.2, y=0.5, labels = A4path1w6)

A1path2w6 <- round(sum(as.numeric((path2w6 <= nw) & (path2w6 >= tff$W)))/(1087/100), digits = 1)

A2path2w6 <- round(sum(as.numeric((path2w6 >= nw)&(path2w6 >= tff$W)))/(1087/100), digits = 1)

A3path2w6 <- round(sum(as.numeric((path2w6 >= nw)&(path2w6 <= tff$W)))/(1087/100), digits=1)

A4path2w6 <- round(sum(as.numeric((path2w6 <= nw)&(path2w6 <= tff$W)))/(1087/100), digits = 1)

plot(nw, path2w6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w6)
text(x=0.8, y=0.5, labels = A2path2w6)
text(x=0.5, y=0.2, labels = A3path2w6)
text(x=0.2, y=0.5, labels = A4path2w6)

A1epsolw6 <- round(sum(as.numeric((epsolw6 <= nw) & (epsolw6 >= tff$W)))/(1087/100), digits = 1)

A2epsolw6 <- round(sum(as.numeric((epsolw6 >= nw)&(epsolw6 >= tff$W)))/(1087/100), digits = 1)

A3epsolw6 <- round(sum(as.numeric((epsolw6 >= nw)&(epsolw6 <= tff$W)))/(1087/100), digits=1)

A4epsolw6 <- round(sum(as.numeric((epsolw6 <= nw)&(epsolw6 <= tff$W)))/(1087/100), digits = 1)


plot(nw, epsolw6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw6)
text(x=0.8, y=0.5, labels = A2epsolw6)
text(x=0.5, y=0.2, labels = A3epsolw6)
text(x=0.2, y=0.5, labels = A4epsolw6)



#####w7  frequency threshold 54, raw consistency threshold 0.75#####


ttw7 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.75, n.cut=54, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw7

#parsimonious solution
psw7 <- eqmcc(ttw7, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw7

#exclude contradictory assumptions & rows: OP*SM*CM + SP*SM*CM + TP*SM*CM (contradict statement of sufficiency)

ttw7new <- ttw7
suf
rows <- findRows("OP*SM*CM + SP*SM*CM + TP*SM*CM", ttw7new, remainders = FALSE)
rows

ttw7new$tt[as.character(rows), "OUT"] <- 0
ttw7new


#conservative solution
csw7 <- eqmcc(ttw7new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw7

#enhanced parsimonious solution (all simplifying assumptions)

epsw7 <- eqmcc(ttw7new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw7


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  sm
epsolw7 <- compute("sm", data=tff)

obsepsolw7 <- as.numeric(epsolw7 <= nw)
zepsolw7 <- round(((((sum(obsepsolw7 == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits=1)
zepsolw7
pzepsolw7 <- 2*pnorm(-abs(zepsolw7))
pzepsolw7


#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw7 <- eqmcc(ttw7new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw7


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: sm
path1w7 <- compute("sm ", data=tff)


A1path1w7 <- round(sum(as.numeric((path1w7 <= nw) & (path1w7 >= tff$W)))/(1087/100), digits = 1)

A2path1w7 <- round(sum(as.numeric((path1w7 >= nw)&(path1w7 >= tff$W)))/(1087/100), digits = 1)

A3path1w7 <- round(sum(as.numeric((path1w7 >= nw)&(path1w7 <= tff$W)))/(1087/100), digits=1)

A4path1w7 <- round(sum(as.numeric((path1w7 <= nw)&(path1w7 <= tff$W)))/(1087/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w7, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w7)
text(x=0.8, y=0.5, labels = A2path1w7)
text(x=0.5, y=0.2, labels = A3path1w7)
text(x=0.2, y=0.5, labels = A4path1w7)





###############################THEORY EVALUATION COVERAGE ############

tff <- read.csv("tff_d2_c1.csv", row.names=1, header = TRUE, sep = ";",  dec = ",")
describe(tff)
nw <- 1-tff$W



## Theory evaluation for hypothesis 2####

## calculate intersections

Tw <- "sp*tp*op"
Tw
tw <- deMorgan("sp*tp*op")
tw
Sw <- "op*sm*CM + sp*sm*CM + sp*tp*op*CM + sp*tp*OP*SM*cm + SP*tp*op*SM*cm"
Sw
sw <- deMorgan("op*sm*CM + sp*sm*CM + sp*tp*op*CM + sp*tp*OP*SM*cm + SP*tp*op*SM*cm")
sw

TSw <- intersection("sp*tp*op", "op*sm*CM + sp*sm*CM + sp*tp*op*CM + sp*tp*OP*SM*cm + SP*tp*op*SM*cm", snames = "SP, TP, OP, SM, CM")
TSw

tw
Sw
tSw <- intersection("OP + SP + TP", "op*sm*CM + sp*sm*CM + sp*tp*op*CM + sp*tp*OP*SM*cm + SP*tp*op*SM*cm", snames = "SP, TP, OP, SM, CM")
tSw

sw
tsw <- intersection("OP + SP + TP", "cm*op*sp + cm*sm + cm*TP + CM*OP*SM + OP*SP + CM*SM*SP + SM*TP", snames = "SP, TP, OP, SM, CM")
tsw

Tw
Tsw <- intersection("sp*tp*op", "cm*op*sp + cm*sm + cm*TP + CM*OP*SM + OP*SP + CM*SM*SP + SM*TP", snames = "SP, TP, OP, SM, CM")
Tsw

#calculate membership in intersections



TSw
TSw1 <- compute("sp*tp*op*sm*CM", data=tff)
TSw2 <- compute("sp*tp*op*CM", data=tff)
TSw <- compute("sp*tp*op*sm*CM + sp*tp*op*CM", data=tff)

Tsw
Tsw1 <- compute("sp*tp*op*cm", data=tff)
Tsw2 <- compute("sp*tp*op*sm*cm", data=tff)
Tsw <- compute("sp*tp*op*cm + sp*tp*op*sm*cm", data=tff)


tSw
tSw1 <- compute("sp*OP*sm*CM", data=tff)
tSw2 <- compute("sp*tp*OP*SM*cm", data=tff)
tSw3 <- compute("SP*op*sm*CM", data=tff)
tSw4 <- compute("SP*tp*op*SM*cm", data=tff)
tSw5 <- compute("TP*op*sm*CM", data=tff)
tSw6 <- compute("sp*TP*sm*CM", data=tff)
tSw <- compute("sp*OP*sm*CM + sp*tp*OP*SM*cm + SP*op*sm*CM + SP*tp*op*SM*cm + TP*op*sm*CM + sp*TP*sm*CM", data=tff)

tsw
tsw1 <- compute("OP*sm*cm", data=tff)
tsw2 <- compute("TP*OP*cm", data=tff)
tsw3 <- compute("OP*SM*CM", data=tff)
tsw4 <- compute("SP*OP", data=tff)
tsw5 <- compute("SP*OP*SM*CM", data=tff)
tsw6 <- compute("TP*OP*SM", data=tff)
tsw7 <- compute("SP*sm*cm", data=tff)
tsw8 <- compute("SP*TP*cm", data=tff)
tsw9 <- compute("SP*SM*CM", data=tff)
tsw10 <- compute("SP*TP*SM", data=tff)
tsw11 <- compute("sp*TP*op*cm", data=tff)
tsw12 <- compute("TP*sm*cm", data=tff)
tsw13 <- compute("TP*cm", data=tff)
tsw14 <- compute("TP*OP*SM*CM", data=tff)
tsw15 <- compute("SP*TP*OP", data=tff)
tsw16 <- compute("SP*TP*SM*CM", data=tff)
tsw17 <- compute("TP*SM", data=tff)
tsw <- compute("OP*sm*cm + TP*OP*cm + OP*SM*CM + SP*OP + SP*OP*SM*CM + TP*OP*SM + SP*sm*cm + SP*TP*cm + SP*SM*CM + SP*TP*SM + sp*TP*op*cm + TP*sm*cm + TP*cm + TP*OP*SM*CM + SP*TP*OP + SP*TP*SM*CM + TP*SM", data=tff)



#check for logical remainders
TSw
sum(as.numeric(TSw1 > 0.5))
sum(as.numeric(TSw2 > 0.5))
sum(as.numeric(TSw > 0.5))

Tsw
sum(as.numeric(Tsw1 > 0.5))
sum(as.numeric(Tsw2 > 0.5))
sum(as.numeric(Tsw > 0.5))

sum(as.numeric(tSw1 > 0.5))
sum(as.numeric(tSw2 > 0.5))
sum(as.numeric(tSw3 > 0.5))
sum(as.numeric(tSw4 > 0.5))
sum(as.numeric(tSw5 > 0.5))
sum(as.numeric(tSw6 > 0.5))
sum(as.numeric(tSw > 0.5))

sum(as.numeric(tsw1 > 0.5))
sum(as.numeric(tsw2 > 0.5))
sum(as.numeric(tsw3 > 0.5))
sum(as.numeric(tsw4 > 0.5))
sum(as.numeric(tsw5 > 0.5))
sum(as.numeric(tsw6 > 0.5))
sum(as.numeric(tsw7 > 0.5))
sum(as.numeric(tsw8 > 0.5))
sum(as.numeric(tsw9 > 0.5))
sum(as.numeric(tsw10 > 0.5))
sum(as.numeric(tsw11 > 0.5))
sum(as.numeric(tsw12 > 0.5))
sum(as.numeric(tsw13 > 0.5))
sum(as.numeric(tsw14 > 0.5))
sum(as.numeric(tsw15 > 0.5))
sum(as.numeric(tsw16 > 0.5))
sum(as.numeric(tsw17 > 0.5))
sum(as.numeric(tsw > 0.5))


#calculate % with W and w

TSwW <- round(sum(as.numeric((TSw > 0.5) & (tff$W > 0.5)))/(1087/100), digits = 1)
TSwW

TSww <- round(sum(as.numeric((TSw > 0.5 )&(tff$W < 0.5)))/(1087/100), digits = 1)
TSww

TswW <- round(sum(as.numeric((Tsw > 0.5) & (tff$W > 0.5)))/(1087/100), digits = 1)
TswW

Tsww <- round(sum(as.numeric((Tsw > 0.5 )&(tff$W < 0.5)))/(1087/100), digits = 1)
Tsww

tswW <- round(sum(as.numeric((tsw > 0.5) & (tff$W > 0.5)))/(1087/100), digits = 1)
tswW

tsww <- round(sum(as.numeric((tsw > 0.5 )&(tff$W < 0.5)))/(1087/100), digits = 1)
tsww

tSwW <- round(sum(as.numeric((tSw > 0.5) & (tff$W > 0.5)))/(1087/100), digits = 1)
tSwW

tSww <- round(sum(as.numeric((tSw > 0.5 )&(tff$W < 0.5)))/(1087/100), digits = 1)
tSww

#check for mistakes
fulltew <- sum(TSwW, TSww, tSwW, tSww, TswW, Tsww, tswW, tsww)
fulltew


##########hier stehengeblieben######
####Theory evaluation for hypothesis 3#####
TW <- " OP*SM + OP*CM + SP*SM + SP*CM + TP*SM + TP*CM "
TW
tW <- deMorgan("OP*SM + OP*CM + SP*SM + SP*CM + TP*SM + TP*CM ")
tW
SW <- "OP*SM*CM + SP*SM*CM + TP*SM*CM"
SW
sW <- deMorgan("OP*SM*CM + SP*SM*CM + TP*SM*CM")
sW

TSW <- intersection("OP*SM + OP*CM + SP*SM + SP*CM + TP*SM + TP*CM ", " OP*SM*CM + SP*SM*CM + TP*SM*CM ", snames = "SP, TP, OP, SM, CM")
TSW

tW
SW
tSW <- intersection("op*sp*tp + cm*sm", "OP*SM*CM + SP*SM*CM + TP*SM*CM", snames = "SP, TP, OP, SM, CM")
tSW


sW
tW
tsW <- intersection("op*sp*tp + cm*sm", "cm + op*sp*tp + sm", snames = "SP, TP, OP, SM, CM")
tsW

TW
sW
TsW <- intersection(" OP*SM + OP*CM + SP*SM + SP*CM + TP*SM + TP*CM ", "cm + op*sp*tp + sm", snames = "SP, TP, OP, SM, CM")
TsW


#calculate membership in intersections

TSW
TSW1 <- compute("OP*SM*CM", data=tff)
TSW2 <- compute("SP*OP*SM*CM", data=tff)
TSW3 <- compute("TP*OP*SM*CM", data=tff)
TSW4 <- compute("SP*SM*CM", data=tff)
TSW5 <- compute("SP*TP*SM*CM", data=tff)
TSW6 <- compute("TP*SM*CM", data=tff)

TSW <- compute("OP*SM*CM + SP*OP*SM*CM + TP*OP*SM*CM + SP*SM*CM + SP*TP*SM*CM + TP*SM*CM", data=tff)

TsW
TsW1 <- compute("OP*SM*cm", data=tff)
TsW2 <- compute("OP*sm*CM", data=tff)
TsW3 <- compute("SP*SM*cm", data=tff)
TsW4 <- compute("SP*sm*CM", data=tff)
TsW5 <- compute("TP*SM*cm", data=tff)
TsW6 <- compute("TP*sm*CM", data=tff)

TsW <- compute("OP*SM*cm + OP*sm*CM + SP*SM*cm + SP*sm*CM + TP*SM*cm + TP*sm*CM", data=tff)

tSW

tsW
tsW1 <- compute("sp*tp*op*cm", data=tff)
tsW2 <- compute("sp*tp*op", data=tff)
tsW3 <- compute("sp*tp*op*sm", data=tff)
tsW4 <- compute("sm*cm", data=tff)
tsW5 <- compute("sp*tp*op*sm*cm", data=tff)

tsW <- compute("sp*tp*op*cm + sp*tp*op + sp*tp*op*sm + sm*cm + sp*tp*op*sm*cm", data=tff)

#check for logical remainders
TSW
sum(as.numeric(TSW1 > 0.5))
sum(as.numeric(TSW2 > 0.5))
sum(as.numeric(TSW3 > 0.5))
sum(as.numeric(TSW4 > 0.5))
sum(as.numeric(TSW5 > 0.5))
sum(as.numeric(TSW6 > 0.5))


sum(as.numeric(TSW > 0.5))

TsW
sum(as.numeric(TsW1 > 0.5))
sum(as.numeric(TsW2 > 0.5))
sum(as.numeric(TsW3 > 0.5))
sum(as.numeric(TsW4 > 0.5))
sum(as.numeric(TsW5 > 0.5))
sum(as.numeric(TsW6 > 0.5))

sum(as.numeric(TsW > 0.5))

tsW
sum(as.numeric(tsW1 > 0.5))
sum(as.numeric(tsW2 > 0.5))
sum(as.numeric(tsW3 > 0.5))
sum(as.numeric(tsW4 > 0.5))
sum(as.numeric(tsW5 > 0.5))

sum(as.numeric(tsW > 0.5))

#calculate % with W and w

TSWW <- round(sum(as.numeric((TSW > 0.5) & (tff$W > 0.5)))/(1087/100), digits = 1)
TSWW

TSWw <- round(sum(as.numeric((TSW > 0.5 )&(tff$W < 0.5)))/(1087/100), digits = 1)
TSWw

TsWW <- round(sum(as.numeric((TsW > 0.5) & (tff$W > 0.5)))/(1087/100), digits = 1)
TsWW

TsWw <- round(sum(as.numeric((TsW > 0.5 )&(tff$W < 0.5)))/(1087/100), digits = 1)
TsWw

tsWW <- round(sum(as.numeric((tsW > 0.5) & (tff$W > 0.5)))/(1087/100), digits = 1)
tsWW

tsWw <- round(sum(as.numeric((tsW > 0.5 )&(tff$W < 0.5)))/(1087/100), digits = 1)
tsWw

tSWW <- round(sum(as.numeric((tSW > 0.5) & (tff$W > 0.5)))/(1087/100), digits = 1)
tSWW

tSWw <- round(sum(as.numeric((tSW > 0.5 )&(tff$W < 0.5)))/(1087/100), digits = 1)
tSWw




#check for mistakes
fulltew <- sum(TSWW, TSWw, tSWW, tSWw, TsWW, TsWw, tsWW, tsWw)
fulltew



